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