 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CGET03( 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 REAL RCOND, RESID 00011 * .. 00012 * .. Array Arguments .. 00013 REAL RWORK( * ) 00014 COMPLEX A( LDA, * ), AINV( LDAINV, * ), 00015 $ WORK( LDWORK, * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * CGET03 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 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 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 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) REAL array, dimension (N) 00049 * 00050 * RCOND (output) REAL 00051 * The reciprocal of the condition number of A, computed as 00052 * ( 1/norm(A) ) / norm(AINV). 00053 * 00054 * RESID (output) REAL 00055 * norm(I - AINV*A) / ( N * norm(A) * norm(AINV) * EPS ) 00056 * 00057 * ===================================================================== 00058 * 00059 * .. Parameters .. 00060 REAL ZERO, ONE 00061 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00062 COMPLEX CZERO, CONE 00063 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00064 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00065 * .. 00066 * .. Local Scalars .. 00067 INTEGER I 00068 REAL AINVNM, ANORM, EPS 00069 * .. 00070 * .. External Functions .. 00071 REAL CLANGE, SLAMCH 00072 EXTERNAL CLANGE, SLAMCH 00073 * .. 00074 * .. External Subroutines .. 00075 EXTERNAL CGEMM 00076 * .. 00077 * .. Intrinsic Functions .. 00078 INTRINSIC REAL 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 = SLAMCH( 'Epsilon' ) 00093 ANORM = CLANGE( '1', N, N, A, LDA, RWORK ) 00094 AINVNM = CLANGE( '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 CGEMM( 'No transpose', 'No transpose', N, N, N, -CONE, 00105 $ AINV, 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 = CLANGE( '1', N, N, WORK, LDWORK, RWORK ) 00113 * 00114 RESID = ( ( RESID*RCOND )/EPS ) / REAL( N ) 00115 * 00116 RETURN 00117 * 00118 * End of CGET03 00119 * 00120 END
1.7.3 
