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