1*> \brief \b ZHET01_AA 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8* Definition: 9* =========== 10* 11* SUBROUTINE ZHET01_AA( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, 12* C, LDC, RWORK, RESID ) 13* 14* .. Scalar Arguments .. 15* CHARACTER UPLO 16* INTEGER LDA, LDAFAC, LDC, N 17* DOUBLE PRECISION RESID 18* .. 19* .. Array Arguments .. 20* INTEGER IPIV( * ) 21* DOUBLE PRECISION RWORK( * ) 22* COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ) 23* .. 24* 25* 26*> \par Purpose: 27* ============= 28*> 29*> \verbatim 30*> 31*> ZHET01_AA reconstructs a hermitian indefinite matrix A from its 32*> block L*D*L' or U*D*U' factorization and computes the residual 33*> norm( C - A ) / ( N * norm(A) * EPS ), 34*> where C is the reconstructed matrix and EPS is the machine epsilon. 35*> \endverbatim 36* 37* Arguments: 38* ========== 39* 40*> \param[in] UPLO 41*> \verbatim 42*> UPLO is CHARACTER*1 43*> Specifies whether the upper or lower triangular part of the 44*> hermitian matrix A is stored: 45*> = 'U': Upper triangular 46*> = 'L': Lower triangular 47*> \endverbatim 48*> 49*> \param[in] N 50*> \verbatim 51*> N is INTEGER 52*> The number of rows and columns of the matrix A. N >= 0. 53*> \endverbatim 54*> 55*> \param[in] A 56*> \verbatim 57*> A is COMPLEX*16 array, dimension (LDA,N) 58*> The original hermitian matrix A. 59*> \endverbatim 60*> 61*> \param[in] LDA 62*> \verbatim 63*> LDA is INTEGER 64*> The leading dimension of the array A. LDA >= max(1,N) 65*> \endverbatim 66*> 67*> \param[in] AFAC 68*> \verbatim 69*> AFAC is COMPLEX*16 array, dimension (LDAFAC,N) 70*> The factored form of the matrix A. AFAC contains the block 71*> diagonal matrix D and the multipliers used to obtain the 72*> factor L or U from the block L*D*L' or U*D*U' factorization 73*> as computed by ZHETRF. 74*> \endverbatim 75*> 76*> \param[in] LDAFAC 77*> \verbatim 78*> LDAFAC is INTEGER 79*> The leading dimension of the array AFAC. LDAFAC >= max(1,N). 80*> \endverbatim 81*> 82*> \param[in] IPIV 83*> \verbatim 84*> IPIV is INTEGER array, dimension (N) 85*> The pivot indices from ZHETRF. 86*> \endverbatim 87*> 88*> \param[out] C 89*> \verbatim 90*> C is COMPLEX*16 array, dimension (LDC,N) 91*> \endverbatim 92*> 93*> \param[in] LDC 94*> \verbatim 95*> LDC is INTEGER 96*> The leading dimension of the array C. LDC >= max(1,N). 97*> \endverbatim 98*> 99*> \param[out] RWORK 100*> \verbatim 101*> RWORK is COMPLEX*16 array, dimension (N) 102*> \endverbatim 103*> 104*> \param[out] RESID 105*> \verbatim 106*> RESID is COMPLEX*16 107*> If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS ) 108*> If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS ) 109*> \endverbatim 110* 111* Authors: 112* ======== 113* 114*> \author Univ. of Tennessee 115*> \author Univ. of California Berkeley 116*> \author Univ. of Colorado Denver 117*> \author NAG Ltd. 118* 119*> \date December 2016 120* 121* 122*> \ingroup complex16_lin 123* 124* ===================================================================== 125 SUBROUTINE ZHET01_AA( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, 126 $ LDC, RWORK, RESID ) 127* 128* -- LAPACK test routine (version 3.7.0) -- 129* -- LAPACK is a software package provided by Univ. of Tennessee, -- 130* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 131* December 2016 132* 133* .. Scalar Arguments .. 134 CHARACTER UPLO 135 INTEGER LDA, LDAFAC, LDC, N 136 DOUBLE PRECISION RESID 137* .. 138* .. Array Arguments .. 139 INTEGER IPIV( * ) 140 DOUBLE PRECISION RWORK( * ) 141 COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ) 142* .. 143* 144* ===================================================================== 145* 146* .. Parameters .. 147 COMPLEX*16 CZERO, CONE 148 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 149 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 150 DOUBLE PRECISION ZERO, ONE 151 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 152* .. 153* .. Local Scalars .. 154 INTEGER I, J 155 DOUBLE PRECISION ANORM, EPS 156* .. 157* .. External Functions .. 158 LOGICAL LSAME 159 DOUBLE PRECISION DLAMCH, ZLANHE 160 EXTERNAL LSAME, DLAMCH, ZLANHE 161* .. 162* .. External Subroutines .. 163 EXTERNAL ZLASET, ZLAVHE 164* .. 165* .. Intrinsic Functions .. 166 INTRINSIC DBLE 167* .. 168* .. Executable Statements .. 169* 170* Quick exit if N = 0. 171* 172 IF( N.LE.0 ) THEN 173 RESID = ZERO 174 RETURN 175 END IF 176* 177* Determine EPS and the norm of A. 178* 179 EPS = DLAMCH( 'Epsilon' ) 180 ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK ) 181* 182* Initialize C to the tridiagonal matrix T. 183* 184 CALL ZLASET( 'Full', N, N, CZERO, CZERO, C, LDC ) 185 CALL ZLACPY( 'F', 1, N, AFAC( 1, 1 ), LDAFAC+1, C( 1, 1 ), LDC+1 ) 186 IF( N.GT.1 ) THEN 187 IF( LSAME( UPLO, 'U' ) ) THEN 188 CALL ZLACPY( 'F', 1, N-1, AFAC( 1, 2 ), LDAFAC+1, C( 1, 2 ), 189 $ LDC+1 ) 190 CALL ZLACPY( 'F', 1, N-1, AFAC( 1, 2 ), LDAFAC+1, C( 2, 1 ), 191 $ LDC+1 ) 192 CALL ZLACGV( N-1, C( 2, 1 ), LDC+1 ) 193 ELSE 194 CALL ZLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 1, 2 ), 195 $ LDC+1 ) 196 CALL ZLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 2, 1 ), 197 $ LDC+1 ) 198 CALL ZLACGV( N-1, C( 1, 2 ), LDC+1 ) 199 ENDIF 200* 201* Call ZTRMM to form the product U' * D (or L * D ). 202* 203 IF( LSAME( UPLO, 'U' ) ) THEN 204 CALL ZTRMM( 'Left', UPLO, 'Conjugate transpose', 'Unit', 205 $ N-1, N, CONE, AFAC( 1, 2 ), LDAFAC, C( 2, 1 ), 206 $ LDC ) 207 ELSE 208 CALL ZTRMM( 'Left', UPLO, 'No transpose', 'Unit', N-1, N, 209 $ CONE, AFAC( 2, 1 ), LDAFAC, C( 2, 1 ), LDC ) 210 END IF 211* 212* Call ZTRMM again to multiply by U (or L ). 213* 214 IF( LSAME( UPLO, 'U' ) ) THEN 215 CALL ZTRMM( 'Right', UPLO, 'No transpose', 'Unit', N, N-1, 216 $ CONE, AFAC( 1, 2 ), LDAFAC, C( 1, 2 ), LDC ) 217 ELSE 218 CALL ZTRMM( 'Right', UPLO, 'Conjugate transpose', 'Unit', N, 219 $ N-1, CONE, AFAC( 2, 1 ), LDAFAC, C( 1, 2 ), 220 $ LDC ) 221 END IF 222* 223* Apply hermitian pivots 224* 225 DO J = N, 1, -1 226 I = IPIV( J ) 227 IF( I.NE.J ) 228 $ CALL ZSWAP( N, C( J, 1 ), LDC, C( I, 1 ), LDC ) 229 END DO 230 DO J = N, 1, -1 231 I = IPIV( J ) 232 IF( I.NE.J ) 233 $ CALL ZSWAP( N, C( 1, J ), 1, C( 1, I ), 1 ) 234 END DO 235 ENDIF 236* 237* 238* Compute the difference C - A . 239* 240 IF( LSAME( UPLO, 'U' ) ) THEN 241 DO J = 1, N 242 DO I = 1, J 243 C( I, J ) = C( I, J ) - A( I, J ) 244 END DO 245 END DO 246 ELSE 247 DO J = 1, N 248 DO I = J, N 249 C( I, J ) = C( I, J ) - A( I, J ) 250 END DO 251 END DO 252 END IF 253* 254* Compute norm( C - A ) / ( N * norm(A) * EPS ) 255* 256 RESID = ZLANHE( '1', UPLO, N, C, LDC, RWORK ) 257* 258 IF( ANORM.LE.ZERO ) THEN 259 IF( RESID.NE.ZERO ) 260 $ RESID = ONE / EPS 261 ELSE 262 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS 263 END IF 264* 265 RETURN 266* 267* End of ZHET01_AA 268* 269 END 270