00001 SUBROUTINE CSYT03( 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 RWORK( * ) 00015 COMPLEX A( LDA, * ), AINV( LDAINV, * ), 00016 $ WORK( LDWORK, * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * CSYT03 computes the residual for a complex symmetric matrix times 00023 * its 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 * complex symmetric 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 array, dimension (LDA,N) 00040 * The original complex symmetric 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 array, dimension (LDAINV,N) 00046 * On entry, the inverse of the matrix A, stored as a symmetric 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 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) REAL array, dimension (N) 00062 * 00063 * RCOND (output) REAL 00064 * The reciprocal of the condition number of A, computed as 00065 * RCOND = 1/ (norm(A) * norm(AINV)). 00066 * 00067 * RESID (output) REAL 00068 * norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS ) 00069 * 00070 * ===================================================================== 00071 * 00072 * 00073 * .. Parameters .. 00074 REAL ZERO, ONE 00075 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00076 COMPLEX CZERO, CONE 00077 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00078 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00079 * .. 00080 * .. Local Scalars .. 00081 INTEGER I, J 00082 REAL AINVNM, ANORM, EPS 00083 * .. 00084 * .. External Functions .. 00085 LOGICAL LSAME 00086 REAL CLANGE, CLANSY, SLAMCH 00087 EXTERNAL LSAME, CLANGE, CLANSY, SLAMCH 00088 * .. 00089 * .. External Subroutines .. 00090 EXTERNAL CSYMM 00091 * .. 00092 * .. Intrinsic Functions .. 00093 INTRINSIC REAL 00094 * .. 00095 * .. Executable Statements .. 00096 * 00097 * Quick exit if N = 0 00098 * 00099 IF( N.LE.0 ) THEN 00100 RCOND = ONE 00101 RESID = ZERO 00102 RETURN 00103 END IF 00104 * 00105 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0. 00106 * 00107 EPS = SLAMCH( 'Epsilon' ) 00108 ANORM = CLANSY( '1', UPLO, N, A, LDA, RWORK ) 00109 AINVNM = CLANSY( '1', UPLO, N, AINV, LDAINV, RWORK ) 00110 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 00111 RCOND = ZERO 00112 RESID = ONE / EPS 00113 RETURN 00114 END IF 00115 RCOND = ( ONE/ANORM ) / AINVNM 00116 * 00117 * Expand AINV into a full matrix and call CSYMM to multiply 00118 * AINV on the left by A (store the result in WORK). 00119 * 00120 IF( LSAME( UPLO, 'U' ) ) THEN 00121 DO 20 J = 1, N 00122 DO 10 I = 1, J - 1 00123 AINV( J, I ) = AINV( I, J ) 00124 10 CONTINUE 00125 20 CONTINUE 00126 ELSE 00127 DO 40 J = 1, N 00128 DO 30 I = J + 1, N 00129 AINV( J, I ) = AINV( I, J ) 00130 30 CONTINUE 00131 40 CONTINUE 00132 END IF 00133 CALL CSYMM( 'Left', UPLO, N, N, -CONE, A, LDA, AINV, LDAINV, 00134 $ CZERO, WORK, LDWORK ) 00135 * 00136 * Add the identity matrix to WORK . 00137 * 00138 DO 50 I = 1, N 00139 WORK( I, I ) = WORK( I, I ) + CONE 00140 50 CONTINUE 00141 * 00142 * Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS) 00143 * 00144 RESID = CLANGE( '1', N, N, WORK, LDWORK, RWORK ) 00145 * 00146 RESID = ( ( RESID*RCOND )/EPS ) / REAL( N ) 00147 * 00148 RETURN 00149 * 00150 * End of CSYT03 00151 * 00152 END
1.7.2