1*> \brief \b ZPOT01 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 ZPOT01( UPLO, N, A, LDA, AFAC, LDAFAC, RWORK, RESID ) 12* 13* .. Scalar Arguments .. 14* CHARACTER UPLO 15* INTEGER LDA, LDAFAC, N 16* DOUBLE PRECISION RESID 17* .. 18* .. Array Arguments .. 19* DOUBLE PRECISION RWORK( * ) 20* COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ) 21* .. 22* 23* 24*> \par Purpose: 25* ============= 26*> 27*> \verbatim 28*> 29*> ZPOT01 reconstructs a Hermitian positive definite matrix A from 30*> its L*L' or U'*U factorization and computes the residual 31*> norm( L*L' - A ) / ( N * norm(A) * EPS ) or 32*> norm( U'*U - A ) / ( N * norm(A) * EPS ), 33*> where EPS is the machine epsilon, L' is the conjugate transpose of L, 34*> and U' is the conjugate transpose of U. 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,out] AFAC 68*> \verbatim 69*> AFAC is COMPLEX*16 array, dimension (LDAFAC,N) 70*> On entry, the factor L or U from the L * L**H or U**H * U 71*> factorization of A. 72*> Overwritten with the reconstructed matrix, and then with 73*> the difference L * L**H - A (or U**H * U - A). 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[out] RWORK 83*> \verbatim 84*> RWORK is DOUBLE PRECISION array, dimension (N) 85*> \endverbatim 86*> 87*> \param[out] RESID 88*> \verbatim 89*> RESID is DOUBLE PRECISION 90*> If UPLO = 'L', norm(L * L**H - A) / ( N * norm(A) * EPS ) 91*> If UPLO = 'U', norm(U**H * U - A) / ( N * norm(A) * EPS ) 92*> \endverbatim 93* 94* Authors: 95* ======== 96* 97*> \author Univ. of Tennessee 98*> \author Univ. of California Berkeley 99*> \author Univ. of Colorado Denver 100*> \author NAG Ltd. 101* 102*> \ingroup complex16_lin 103* 104* ===================================================================== 105 SUBROUTINE ZPOT01( UPLO, N, A, LDA, AFAC, LDAFAC, RWORK, RESID ) 106* 107* -- LAPACK test routine -- 108* -- LAPACK is a software package provided by Univ. of Tennessee, -- 109* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 110* 111* .. Scalar Arguments .. 112 CHARACTER UPLO 113 INTEGER LDA, LDAFAC, N 114 DOUBLE PRECISION RESID 115* .. 116* .. Array Arguments .. 117 DOUBLE PRECISION RWORK( * ) 118 COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ) 119* .. 120* 121* ===================================================================== 122* 123* .. Parameters .. 124 DOUBLE PRECISION ZERO, ONE 125 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 126* .. 127* .. Local Scalars .. 128 INTEGER I, J, K 129 DOUBLE PRECISION ANORM, EPS, TR 130 COMPLEX*16 TC 131* .. 132* .. External Functions .. 133 LOGICAL LSAME 134 DOUBLE PRECISION DLAMCH, ZLANHE 135 COMPLEX*16 ZDOTC 136 EXTERNAL LSAME, DLAMCH, ZLANHE, ZDOTC 137* .. 138* .. External Subroutines .. 139 EXTERNAL ZHER, ZSCAL, ZTRMV 140* .. 141* .. Intrinsic Functions .. 142 INTRINSIC DBLE, DIMAG 143* .. 144* .. Executable Statements .. 145* 146* Quick exit if N = 0. 147* 148 IF( N.LE.0 ) THEN 149 RESID = ZERO 150 RETURN 151 END IF 152* 153* Exit with RESID = 1/EPS if ANORM = 0. 154* 155 EPS = DLAMCH( 'Epsilon' ) 156 ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK ) 157 IF( ANORM.LE.ZERO ) THEN 158 RESID = ONE / EPS 159 RETURN 160 END IF 161* 162* Check the imaginary parts of the diagonal elements and return with 163* an error code if any are nonzero. 164* 165 DO 10 J = 1, N 166 IF( DIMAG( AFAC( J, J ) ).NE.ZERO ) THEN 167 RESID = ONE / EPS 168 RETURN 169 END IF 170 10 CONTINUE 171* 172* Compute the product U**H * U, overwriting U. 173* 174 IF( LSAME( UPLO, 'U' ) ) THEN 175 DO 20 K = N, 1, -1 176* 177* Compute the (K,K) element of the result. 178* 179 TR = ZDOTC( K, AFAC( 1, K ), 1, AFAC( 1, K ), 1 ) 180 AFAC( K, K ) = TR 181* 182* Compute the rest of column K. 183* 184 CALL ZTRMV( 'Upper', 'Conjugate', 'Non-unit', K-1, AFAC, 185 $ LDAFAC, AFAC( 1, K ), 1 ) 186* 187 20 CONTINUE 188* 189* Compute the product L * L**H, overwriting L. 190* 191 ELSE 192 DO 30 K = N, 1, -1 193* 194* Add a multiple of column K of the factor L to each of 195* columns K+1 through N. 196* 197 IF( K+1.LE.N ) 198 $ CALL ZHER( 'Lower', N-K, ONE, AFAC( K+1, K ), 1, 199 $ AFAC( K+1, K+1 ), LDAFAC ) 200* 201* Scale column K by the diagonal element. 202* 203 TC = AFAC( K, K ) 204 CALL ZSCAL( N-K+1, TC, AFAC( K, K ), 1 ) 205* 206 30 CONTINUE 207 END IF 208* 209* Compute the difference L * L**H - A (or U**H * U - A). 210* 211 IF( LSAME( UPLO, 'U' ) ) THEN 212 DO 50 J = 1, N 213 DO 40 I = 1, J - 1 214 AFAC( I, J ) = AFAC( I, J ) - A( I, J ) 215 40 CONTINUE 216 AFAC( J, J ) = AFAC( J, J ) - DBLE( A( J, J ) ) 217 50 CONTINUE 218 ELSE 219 DO 70 J = 1, N 220 AFAC( J, J ) = AFAC( J, J ) - DBLE( A( J, J ) ) 221 DO 60 I = J + 1, N 222 AFAC( I, J ) = AFAC( I, J ) - A( I, J ) 223 60 CONTINUE 224 70 CONTINUE 225 END IF 226* 227* Compute norm(L*U - A) / ( N * norm(A) * EPS ) 228* 229 RESID = ZLANHE( '1', UPLO, N, AFAC, LDAFAC, RWORK ) 230* 231 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS 232* 233 RETURN 234* 235* End of ZPOT01 236* 237 END 238