 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE ZGET03( N, A, LDA, AINV, LDAINV, WORK, LDWORK, RWORK, 00002 $ 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 INTEGER LDA, LDAINV, LDWORK, N 00010 DOUBLE PRECISION RCOND, RESID 00011 * .. 00012 * .. Array Arguments .. 00013 DOUBLE PRECISION RWORK( * ) 00014 COMPLEX*16 A( LDA, * ), AINV( LDAINV, * ), 00015 $ WORK( LDWORK, * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * ZGET03 computes the residual for a general matrix times its inverse: 00022 * norm( I - AINV*A ) / ( N * norm(A) * norm(AINV) * EPS ), 00023 * where EPS is the machine epsilon. 00024 * 00025 * Arguments 00026 * ========== 00027 * 00028 * N (input) INTEGER 00029 * The number of rows and columns of the matrix A. N >= 0. 00030 * 00031 * A (input) COMPLEX*16 array, dimension (LDA,N) 00032 * The original N x N matrix A. 00033 * 00034 * LDA (input) INTEGER 00035 * The leading dimension of the array A. LDA >= max(1,N). 00036 * 00037 * AINV (input) COMPLEX*16 array, dimension (LDAINV,N) 00038 * The inverse of the matrix A. 00039 * 00040 * LDAINV (input) INTEGER 00041 * The leading dimension of the array AINV. LDAINV >= max(1,N). 00042 * 00043 * WORK (workspace) COMPLEX*16 array, dimension (LDWORK,N) 00044 * 00045 * LDWORK (input) INTEGER 00046 * The leading dimension of the array WORK. LDWORK >= max(1,N). 00047 * 00048 * RWORK (workspace) DOUBLE PRECISION array, dimension (N) 00049 * 00050 * RCOND (output) DOUBLE PRECISION 00051 * The reciprocal of the condition number of A, computed as 00052 * ( 1/norm(A) ) / norm(AINV). 00053 * 00054 * RESID (output) DOUBLE PRECISION 00055 * norm(I - AINV*A) / ( N * norm(A) * norm(AINV) * EPS ) 00056 * 00057 * ===================================================================== 00058 * 00059 * .. Parameters .. 00060 DOUBLE PRECISION ZERO, ONE 00061 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00062 COMPLEX*16 CZERO, CONE 00063 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00064 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 00065 * .. 00066 * .. Local Scalars .. 00067 INTEGER I 00068 DOUBLE PRECISION AINVNM, ANORM, EPS 00069 * .. 00070 * .. External Functions .. 00071 DOUBLE PRECISION DLAMCH, ZLANGE 00072 EXTERNAL DLAMCH, ZLANGE 00073 * .. 00074 * .. External Subroutines .. 00075 EXTERNAL ZGEMM 00076 * .. 00077 * .. Intrinsic Functions .. 00078 INTRINSIC DBLE 00079 * .. 00080 * .. Executable Statements .. 00081 * 00082 * Quick exit if N = 0. 00083 * 00084 IF( N.LE.0 ) THEN 00085 RCOND = ONE 00086 RESID = ZERO 00087 RETURN 00088 END IF 00089 * 00090 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0. 00091 * 00092 EPS = DLAMCH( 'Epsilon' ) 00093 ANORM = ZLANGE( '1', N, N, A, LDA, RWORK ) 00094 AINVNM = ZLANGE( '1', N, N, AINV, LDAINV, RWORK ) 00095 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 00096 RCOND = ZERO 00097 RESID = ONE / EPS 00098 RETURN 00099 END IF 00100 RCOND = ( ONE / ANORM ) / AINVNM 00101 * 00102 * Compute I - A * AINV 00103 * 00104 CALL ZGEMM( 'No transpose', 'No transpose', N, N, N, -CONE, AINV, 00105 $ LDAINV, A, LDA, CZERO, WORK, LDWORK ) 00106 DO 10 I = 1, N 00107 WORK( I, I ) = CONE + WORK( I, I ) 00108 10 CONTINUE 00109 * 00110 * Compute norm(I - AINV*A) / (N * norm(A) * norm(AINV) * EPS) 00111 * 00112 RESID = ZLANGE( '1', N, N, WORK, LDWORK, RWORK ) 00113 * 00114 RESID = ( ( RESID*RCOND ) / EPS ) / DBLE( N ) 00115 * 00116 RETURN 00117 * 00118 * End of ZGET03 00119 * 00120 END
1.7.3 
