1*> \brief \b ZHET01_3 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_3( UPLO, N, A, LDA, AFAC, LDAFAC, E, IPIV, C, 12* 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* E( * ) 24* .. 25* 26* 27*> \par Purpose: 28* ============= 29*> 30*> \verbatim 31*> 32*> ZHET01_3 reconstructs a Hermitian indefinite matrix A from its 33*> block L*D*L' or U*D*U' factorization computed by ZHETRF_RK 34*> (or ZHETRF_BK) and computes the residual 35*> norm( C - A ) / ( N * norm(A) * EPS ), 36*> where C is the reconstructed matrix and EPS is the machine epsilon. 37*> \endverbatim 38* 39* Arguments: 40* ========== 41* 42*> \param[in] UPLO 43*> \verbatim 44*> UPLO is CHARACTER*1 45*> Specifies whether the upper or lower triangular part of the 46*> Hermitian matrix A is stored: 47*> = 'U': Upper triangular 48*> = 'L': Lower triangular 49*> \endverbatim 50*> 51*> \param[in] N 52*> \verbatim 53*> N is INTEGER 54*> The number of rows and columns of the matrix A. N >= 0. 55*> \endverbatim 56*> 57*> \param[in] A 58*> \verbatim 59*> A is COMPLEX*16 array, dimension (LDA,N) 60*> The original Hermitian matrix A. 61*> \endverbatim 62*> 63*> \param[in] LDA 64*> \verbatim 65*> LDA is INTEGER 66*> The leading dimension of the array A. LDA >= max(1,N) 67*> \endverbatim 68*> 69*> \param[in] AFAC 70*> \verbatim 71*> AFAC is COMPLEX*16 array, dimension (LDAFAC,N) 72*> Diagonal of the block diagonal matrix D and factors U or L 73*> as computed by ZHETRF_RK and ZHETRF_BK: 74*> a) ONLY diagonal elements of the Hermitian block diagonal 75*> matrix D on the diagonal of A, i.e. D(k,k) = A(k,k); 76*> (superdiagonal (or subdiagonal) elements of D 77*> should be provided on entry in array E), and 78*> b) If UPLO = 'U': factor U in the superdiagonal part of A. 79*> If UPLO = 'L': factor L in the subdiagonal part of A. 80*> \endverbatim 81*> 82*> \param[in] LDAFAC 83*> \verbatim 84*> LDAFAC is INTEGER 85*> The leading dimension of the array AFAC. 86*> LDAFAC >= max(1,N). 87*> \endverbatim 88*> 89*> \param[in] E 90*> \verbatim 91*> E is COMPLEX*16 array, dimension (N) 92*> On entry, contains the superdiagonal (or subdiagonal) 93*> elements of the Hermitian block diagonal matrix D 94*> with 1-by-1 or 2-by-2 diagonal blocks, where 95*> If UPLO = 'U': E(i) = D(i-1,i),i=2:N, E(1) not referenced; 96*> If UPLO = 'L': E(i) = D(i+1,i),i=1:N-1, E(N) not referenced. 97*> \endverbatim 98*> 99*> \param[in] IPIV 100*> \verbatim 101*> IPIV is INTEGER array, dimension (N) 102*> The pivot indices from ZHETRF_RK (or ZHETRF_BK). 103*> \endverbatim 104*> 105*> \param[out] C 106*> \verbatim 107*> C is COMPLEX*16 array, dimension (LDC,N) 108*> \endverbatim 109*> 110*> \param[in] LDC 111*> \verbatim 112*> LDC is INTEGER 113*> The leading dimension of the array C. LDC >= max(1,N). 114*> \endverbatim 115*> 116*> \param[out] RWORK 117*> \verbatim 118*> RWORK is DOUBLE PRECISION array, dimension (N) 119*> \endverbatim 120*> 121*> \param[out] RESID 122*> \verbatim 123*> RESID is DOUBLE PRECISION 124*> If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS ) 125*> If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS ) 126*> \endverbatim 127* 128* Authors: 129* ======== 130* 131*> \author Univ. of Tennessee 132*> \author Univ. of California Berkeley 133*> \author Univ. of Colorado Denver 134*> \author NAG Ltd. 135* 136*> \ingroup complex16_lin 137* 138* ===================================================================== 139 SUBROUTINE ZHET01_3( UPLO, N, A, LDA, AFAC, LDAFAC, E, IPIV, C, 140 $ LDC, RWORK, RESID ) 141* 142* -- LAPACK test routine -- 143* -- LAPACK is a software package provided by Univ. of Tennessee, -- 144* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 145* 146* .. Scalar Arguments .. 147 CHARACTER UPLO 148 INTEGER LDA, LDAFAC, LDC, N 149 DOUBLE PRECISION RESID 150* .. 151* .. Array Arguments .. 152 INTEGER IPIV( * ) 153 DOUBLE PRECISION RWORK( * ) 154 COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ), 155 $ E( * ) 156* .. 157* 158* ===================================================================== 159* 160* .. Parameters .. 161 DOUBLE PRECISION ZERO, ONE 162 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 163 COMPLEX*16 CZERO, CONE 164 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 165 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 166* .. 167* .. Local Scalars .. 168 INTEGER I, INFO, J 169 DOUBLE PRECISION ANORM, EPS 170* .. 171* .. External Functions .. 172 LOGICAL LSAME 173 DOUBLE PRECISION ZLANHE, DLAMCH 174 EXTERNAL LSAME, ZLANHE, DLAMCH 175* .. 176* .. External Subroutines .. 177 EXTERNAL ZLASET, ZLAVHE_ROOK, ZSYCONVF_ROOK 178* .. 179* .. Intrinsic Functions .. 180 INTRINSIC DIMAG, DBLE 181* .. 182* .. Executable Statements .. 183* 184* Quick exit if N = 0. 185* 186 IF( N.LE.0 ) THEN 187 RESID = ZERO 188 RETURN 189 END IF 190* 191* a) Revert to multiplyers of L 192* 193 CALL ZSYCONVF_ROOK( UPLO, 'R', N, AFAC, LDAFAC, E, IPIV, INFO ) 194* 195* 1) Determine EPS and the norm of A. 196* 197 EPS = DLAMCH( 'Epsilon' ) 198 ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK ) 199* 200* Check the imaginary parts of the diagonal elements and return with 201* an error code if any are nonzero. 202* 203 DO J = 1, N 204 IF( DIMAG( AFAC( J, J ) ).NE.ZERO ) THEN 205 RESID = ONE / EPS 206 RETURN 207 END IF 208 END DO 209* 210* 2) Initialize C to the identity matrix. 211* 212 CALL ZLASET( 'Full', N, N, CZERO, CONE, C, LDC ) 213* 214* 3) Call ZLAVHE_ROOK to form the product D * U' (or D * L' ). 215* 216 CALL ZLAVHE_ROOK( UPLO, 'Conjugate', 'Non-unit', N, N, AFAC, 217 $ LDAFAC, IPIV, C, LDC, INFO ) 218* 219* 4) Call ZLAVHE_RK again to multiply by U (or L ). 220* 221 CALL ZLAVHE_ROOK( UPLO, 'No transpose', 'Unit', N, N, AFAC, 222 $ LDAFAC, IPIV, C, LDC, INFO ) 223* 224* 5) Compute the difference C - A . 225* 226 IF( LSAME( UPLO, 'U' ) ) THEN 227 DO J = 1, N 228 DO I = 1, J - 1 229 C( I, J ) = C( I, J ) - A( I, J ) 230 END DO 231 C( J, J ) = C( J, J ) - DBLE( A( J, J ) ) 232 END DO 233 ELSE 234 DO J = 1, N 235 C( J, J ) = C( J, J ) - DBLE( A( J, J ) ) 236 DO I = J + 1, N 237 C( I, J ) = C( I, J ) - A( I, J ) 238 END DO 239 END DO 240 END IF 241* 242* 6) Compute norm( C - A ) / ( N * norm(A) * EPS ) 243* 244 RESID = ZLANHE( '1', UPLO, N, C, LDC, RWORK ) 245* 246 IF( ANORM.LE.ZERO ) THEN 247 IF( RESID.NE.ZERO ) 248 $ RESID = ONE / EPS 249 ELSE 250 RESID = ( ( RESID/DBLE( N ) )/ANORM ) / EPS 251 END IF 252* 253* b) Convert to factor of L (or U) 254* 255 CALL ZSYCONVF_ROOK( UPLO, 'C', N, AFAC, LDAFAC, E, IPIV, INFO ) 256* 257 RETURN 258* 259* End of ZHET01_3 260* 261 END 262