00001 SUBROUTINE ZPOT03( UPLO, N, A, LDA, AINV, LDAINV, WORK, LDWORK, 00002 $ RWORK, RCOND, RESID ) 00003 * 00004 * -- LAPACK test routine (version 3.1) -- 00005 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 CHARACTER UPLO 00010 INTEGER LDA, LDAINV, LDWORK, N 00011 DOUBLE PRECISION RCOND, RESID 00012 * .. 00013 * .. Array Arguments .. 00014 DOUBLE PRECISION RWORK( * ) 00015 COMPLEX*16 A( LDA, * ), AINV( LDAINV, * ), 00016 $ WORK( LDWORK, * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * ZPOT03 computes the residual for a Hermitian matrix times its 00023 * inverse: 00024 * norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ), 00025 * where EPS is the machine epsilon. 00026 * 00027 * Arguments 00028 * ========== 00029 * 00030 * UPLO (input) CHARACTER*1 00031 * Specifies whether the upper or lower triangular part of the 00032 * Hermitian matrix A is stored: 00033 * = 'U': Upper triangular 00034 * = 'L': Lower triangular 00035 * 00036 * N (input) INTEGER 00037 * The number of rows and columns of the matrix A. N >= 0. 00038 * 00039 * A (input) COMPLEX*16 array, dimension (LDA,N) 00040 * The original Hermitian matrix A. 00041 * 00042 * LDA (input) INTEGER 00043 * The leading dimension of the array A. LDA >= max(1,N) 00044 * 00045 * AINV (input/output) COMPLEX*16 array, dimension (LDAINV,N) 00046 * On entry, the inverse of the matrix A, stored as a Hermitian 00047 * matrix in the same format as A. 00048 * In this version, AINV is expanded into a full matrix and 00049 * multiplied by A, so the opposing triangle of AINV will be 00050 * changed; i.e., if the upper triangular part of AINV is 00051 * stored, the lower triangular part will be used as work space. 00052 * 00053 * LDAINV (input) INTEGER 00054 * The leading dimension of the array AINV. LDAINV >= max(1,N). 00055 * 00056 * WORK (workspace) COMPLEX*16 array, dimension (LDWORK,N) 00057 * 00058 * LDWORK (input) INTEGER 00059 * The leading dimension of the array WORK. LDWORK >= max(1,N). 00060 * 00061 * RWORK (workspace) DOUBLE PRECISION array, dimension (N) 00062 * 00063 * RCOND (output) DOUBLE PRECISION 00064 * The reciprocal of the condition number of A, computed as 00065 * ( 1/norm(A) ) / norm(AINV). 00066 * 00067 * RESID (output) DOUBLE PRECISION 00068 * norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS ) 00069 * 00070 * ===================================================================== 00071 * 00072 * .. Parameters .. 00073 DOUBLE PRECISION ZERO, ONE 00074 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00075 COMPLEX*16 CZERO, CONE 00076 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00077 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 00078 * .. 00079 * .. Local Scalars .. 00080 INTEGER I, J 00081 DOUBLE PRECISION AINVNM, ANORM, EPS 00082 * .. 00083 * .. External Functions .. 00084 LOGICAL LSAME 00085 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHE 00086 EXTERNAL LSAME, DLAMCH, ZLANGE, ZLANHE 00087 * .. 00088 * .. External Subroutines .. 00089 EXTERNAL ZHEMM 00090 * .. 00091 * .. Intrinsic Functions .. 00092 INTRINSIC DBLE, DCONJG 00093 * .. 00094 * .. Executable Statements .. 00095 * 00096 * Quick exit if N = 0. 00097 * 00098 IF( N.LE.0 ) THEN 00099 RCOND = ONE 00100 RESID = ZERO 00101 RETURN 00102 END IF 00103 * 00104 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0. 00105 * 00106 EPS = DLAMCH( 'Epsilon' ) 00107 ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK ) 00108 AINVNM = ZLANHE( '1', UPLO, N, AINV, LDAINV, RWORK ) 00109 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 00110 RCOND = ZERO 00111 RESID = ONE / EPS 00112 RETURN 00113 END IF 00114 RCOND = ( ONE / ANORM ) / AINVNM 00115 * 00116 * Expand AINV into a full matrix and call ZHEMM to multiply 00117 * AINV on the left by A. 00118 * 00119 IF( LSAME( UPLO, 'U' ) ) THEN 00120 DO 20 J = 1, N 00121 DO 10 I = 1, J - 1 00122 AINV( J, I ) = DCONJG( AINV( I, J ) ) 00123 10 CONTINUE 00124 20 CONTINUE 00125 ELSE 00126 DO 40 J = 1, N 00127 DO 30 I = J + 1, N 00128 AINV( J, I ) = DCONJG( AINV( I, J ) ) 00129 30 CONTINUE 00130 40 CONTINUE 00131 END IF 00132 CALL ZHEMM( 'Left', UPLO, N, N, -CONE, A, LDA, AINV, LDAINV, 00133 $ CZERO, WORK, LDWORK ) 00134 * 00135 * Add the identity matrix to WORK . 00136 * 00137 DO 50 I = 1, N 00138 WORK( I, I ) = WORK( I, I ) + CONE 00139 50 CONTINUE 00140 * 00141 * Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS) 00142 * 00143 RESID = ZLANGE( '1', N, N, WORK, LDWORK, RWORK ) 00144 * 00145 RESID = ( ( RESID*RCOND ) / EPS ) / DBLE( N ) 00146 * 00147 RETURN 00148 * 00149 * End of ZPOT03 00150 * 00151 END
1.7.2