1*> \brief \b CHET01_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 CHET01_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* REAL RESID 18* .. 19* .. Array Arguments .. 20* INTEGER IPIV( * ) 21* REAL RWORK( * ) 22* COMPLEX A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ), 23* E( * ) 24* .. 25* 26* 27*> \par Purpose: 28* ============= 29*> 30*> \verbatim 31*> 32*> CHET01_3 reconstructs a Hermitian indefinite matrix A from its 33*> block L*D*L' or U*D*U' factorization computed by CHETRF_RK 34*> (or CHETRF_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 array, dimension (LDAFAC,N) 72*> Diagonal of the block diagonal matrix D and factors U or L 73*> as computed by CHETRF_RK and CHETRF_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 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 CHETRF_RK (or CHETRF_BK). 103*> \endverbatim 104*> 105*> \param[out] C 106*> \verbatim 107*> C is COMPLEX 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 REAL array, dimension (N) 119*> \endverbatim 120*> 121*> \param[out] RESID 122*> \verbatim 123*> RESID is REAL 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*> \date June 2017 137* 138*> \ingroup complex_lin 139* 140* ===================================================================== 141 SUBROUTINE CHET01_3( UPLO, N, A, LDA, AFAC, LDAFAC, E, IPIV, C, 142 $ LDC, RWORK, RESID ) 143* 144* -- LAPACK test routine (version 3.7.1) -- 145* -- LAPACK is a software package provided by Univ. of Tennessee, -- 146* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 147* June 2017 148* 149* .. Scalar Arguments .. 150 CHARACTER UPLO 151 INTEGER LDA, LDAFAC, LDC, N 152 REAL RESID 153* .. 154* .. Array Arguments .. 155 INTEGER IPIV( * ) 156 REAL RWORK( * ) 157 COMPLEX A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ), 158 $ E( * ) 159* .. 160* 161* ===================================================================== 162* 163* .. Parameters .. 164 REAL ZERO, ONE 165 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 166 COMPLEX CZERO, CONE 167 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 168 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 169* .. 170* .. Local Scalars .. 171 INTEGER I, INFO, J 172 REAL ANORM, EPS 173* .. 174* .. External Functions .. 175 LOGICAL LSAME 176 REAL CLANHE, SLAMCH 177 EXTERNAL LSAME, CLANHE, SLAMCH 178* .. 179* .. External Subroutines .. 180 EXTERNAL CLASET, CLAVHE_ROOK, CSYCONVF_ROOK 181* .. 182* .. Intrinsic Functions .. 183 INTRINSIC AIMAG, REAL 184* .. 185* .. Executable Statements .. 186* 187* Quick exit if N = 0. 188* 189 IF( N.LE.0 ) THEN 190 RESID = ZERO 191 RETURN 192 END IF 193* 194* a) Revert to multiplyers of L 195* 196 CALL CSYCONVF_ROOK( UPLO, 'R', N, AFAC, LDAFAC, E, IPIV, INFO ) 197* 198* 1) Determine EPS and the norm of A. 199* 200 EPS = SLAMCH( 'Epsilon' ) 201 ANORM = CLANHE( '1', UPLO, N, A, LDA, RWORK ) 202* 203* Check the imaginary parts of the diagonal elements and return with 204* an error code if any are nonzero. 205* 206 DO J = 1, N 207 IF( AIMAG( AFAC( J, J ) ).NE.ZERO ) THEN 208 RESID = ONE / EPS 209 RETURN 210 END IF 211 END DO 212* 213* 2) Initialize C to the identity matrix. 214* 215 CALL CLASET( 'Full', N, N, CZERO, CONE, C, LDC ) 216* 217* 3) Call CLAVHE_ROOK to form the product D * U' (or D * L' ). 218* 219 CALL CLAVHE_ROOK( UPLO, 'Conjugate', 'Non-unit', N, N, AFAC, 220 $ LDAFAC, IPIV, C, LDC, INFO ) 221* 222* 4) Call ZLAVHE_RK again to multiply by U (or L ). 223* 224 CALL CLAVHE_ROOK( UPLO, 'No transpose', 'Unit', N, N, AFAC, 225 $ LDAFAC, IPIV, C, LDC, INFO ) 226* 227* 5) Compute the difference C - A . 228* 229 IF( LSAME( UPLO, 'U' ) ) THEN 230 DO J = 1, N 231 DO I = 1, J - 1 232 C( I, J ) = C( I, J ) - A( I, J ) 233 END DO 234 C( J, J ) = C( J, J ) - REAL( A( J, J ) ) 235 END DO 236 ELSE 237 DO J = 1, N 238 C( J, J ) = C( J, J ) - REAL( A( J, J ) ) 239 DO I = J + 1, N 240 C( I, J ) = C( I, J ) - A( I, J ) 241 END DO 242 END DO 243 END IF 244* 245* 6) Compute norm( C - A ) / ( N * norm(A) * EPS ) 246* 247 RESID = CLANHE( '1', UPLO, N, C, LDC, RWORK ) 248* 249 IF( ANORM.LE.ZERO ) THEN 250 IF( RESID.NE.ZERO ) 251 $ RESID = ONE / EPS 252 ELSE 253 RESID = ( ( RESID/REAL( N ) )/ANORM ) / EPS 254 END IF 255* 256* b) Convert to factor of L (or U) 257* 258 CALL CSYCONVF_ROOK( UPLO, 'C', N, AFAC, LDAFAC, E, IPIV, INFO ) 259* 260 RETURN 261* 262* End of CHET01_3 263* 264 END 265