 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 REAL FUNCTION SLA_GERCOND ( TRANS, N, A, LDA, AF, LDAF, IPIV, 00002 $ CMODE, C, INFO, WORK, IWORK ) 00003 * 00004 * -- LAPACK routine (version 3.2.1) -- 00005 * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and -- 00006 * -- Jason Riedy of Univ. of California Berkeley. -- 00007 * -- April 2009 -- 00008 * 00009 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00010 * -- Univ. of California Berkeley and NAG Ltd. -- 00011 * 00012 IMPLICIT NONE 00013 * .. 00014 * .. Scalar Arguments .. 00015 CHARACTER TRANS 00016 INTEGER N, LDA, LDAF, INFO, CMODE 00017 * .. 00018 * .. Array Arguments .. 00019 INTEGER IPIV( * ), IWORK( * ) 00020 REAL A( LDA, * ), AF( LDAF, * ), WORK( * ), 00021 $ C( * ) 00022 * .. 00023 * 00024 * Purpose 00025 * ======= 00026 * 00027 * SLA_GERCOND estimates the Skeel condition number of op(A) * op2(C) 00028 * where op2 is determined by CMODE as follows 00029 * CMODE = 1 op2(C) = C 00030 * CMODE = 0 op2(C) = I 00031 * CMODE = -1 op2(C) = inv(C) 00032 * The Skeel condition number cond(A) = norminf( |inv(A)||A| ) 00033 * is computed by computing scaling factors R such that 00034 * diag(R)*A*op2(C) is row equilibrated and computing the standard 00035 * infinity-norm condition number. 00036 * 00037 * Arguments 00038 * ========== 00039 * 00040 * TRANS (input) CHARACTER*1 00041 * Specifies the form of the system of equations: 00042 * = 'N': A * X = B (No transpose) 00043 * = 'T': A**T * X = B (Transpose) 00044 * = 'C': A**H * X = B (Conjugate Transpose = Transpose) 00045 * 00046 * N (input) INTEGER 00047 * The number of linear equations, i.e., the order of the 00048 * matrix A. N >= 0. 00049 * 00050 * A (input) REAL array, dimension (LDA,N) 00051 * On entry, the N-by-N matrix A. 00052 * 00053 * LDA (input) INTEGER 00054 * The leading dimension of the array A. LDA >= max(1,N). 00055 * 00056 * AF (input) REAL array, dimension (LDAF,N) 00057 * The factors L and U from the factorization 00058 * A = P*L*U as computed by SGETRF. 00059 * 00060 * LDAF (input) INTEGER 00061 * The leading dimension of the array AF. LDAF >= max(1,N). 00062 * 00063 * IPIV (input) INTEGER array, dimension (N) 00064 * The pivot indices from the factorization A = P*L*U 00065 * as computed by SGETRF; row i of the matrix was interchanged 00066 * with row IPIV(i). 00067 * 00068 * CMODE (input) INTEGER 00069 * Determines op2(C) in the formula op(A) * op2(C) as follows: 00070 * CMODE = 1 op2(C) = C 00071 * CMODE = 0 op2(C) = I 00072 * CMODE = -1 op2(C) = inv(C) 00073 * 00074 * C (input) REAL array, dimension (N) 00075 * The vector C in the formula op(A) * op2(C). 00076 * 00077 * INFO (output) INTEGER 00078 * = 0: Successful exit. 00079 * i > 0: The ith argument is invalid. 00080 * 00081 * WORK (input) REAL array, dimension (3*N). 00082 * Workspace. 00083 * 00084 * IWORK (input) INTEGER array, dimension (N). 00085 * Workspace.2 00086 * 00087 * ===================================================================== 00088 * 00089 * .. Local Scalars .. 00090 LOGICAL NOTRANS 00091 INTEGER KASE, I, J 00092 REAL AINVNM, TMP 00093 * .. 00094 * .. Local Arrays .. 00095 INTEGER ISAVE( 3 ) 00096 * .. 00097 * .. External Functions .. 00098 LOGICAL LSAME 00099 EXTERNAL LSAME 00100 * .. 00101 * .. External Subroutines .. 00102 EXTERNAL SLACN2, SGETRS, XERBLA 00103 * .. 00104 * .. Intrinsic Functions .. 00105 INTRINSIC ABS, MAX 00106 * .. 00107 * .. Executable Statements .. 00108 * 00109 SLA_GERCOND = 0.0 00110 * 00111 INFO = 0 00112 NOTRANS = LSAME( TRANS, 'N' ) 00113 IF ( .NOT. NOTRANS .AND. .NOT. LSAME(TRANS, 'T') 00114 $ .AND. .NOT. LSAME(TRANS, 'C') ) THEN 00115 INFO = -1 00116 ELSE IF( N.LT.0 ) THEN 00117 INFO = -2 00118 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00119 INFO = -4 00120 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN 00121 INFO = -6 00122 END IF 00123 IF( INFO.NE.0 ) THEN 00124 CALL XERBLA( 'SLA_GERCOND', -INFO ) 00125 RETURN 00126 END IF 00127 IF( N.EQ.0 ) THEN 00128 SLA_GERCOND = 1.0 00129 RETURN 00130 END IF 00131 * 00132 * Compute the equilibration matrix R such that 00133 * inv(R)*A*C has unit 1-norm. 00134 * 00135 IF (NOTRANS) THEN 00136 DO I = 1, N 00137 TMP = 0.0 00138 IF ( CMODE .EQ. 1 ) THEN 00139 DO J = 1, N 00140 TMP = TMP + ABS( A( I, J ) * C( J ) ) 00141 END DO 00142 ELSE IF ( CMODE .EQ. 0 ) THEN 00143 DO J = 1, N 00144 TMP = TMP + ABS( A( I, J ) ) 00145 END DO 00146 ELSE 00147 DO J = 1, N 00148 TMP = TMP + ABS( A( I, J ) / C( J ) ) 00149 END DO 00150 END IF 00151 WORK( 2*N+I ) = TMP 00152 END DO 00153 ELSE 00154 DO I = 1, N 00155 TMP = 0.0 00156 IF ( CMODE .EQ. 1 ) THEN 00157 DO J = 1, N 00158 TMP = TMP + ABS( A( J, I ) * C( J ) ) 00159 END DO 00160 ELSE IF ( CMODE .EQ. 0 ) THEN 00161 DO J = 1, N 00162 TMP = TMP + ABS( A( J, I ) ) 00163 END DO 00164 ELSE 00165 DO J = 1, N 00166 TMP = TMP + ABS( A( J, I ) / C( J ) ) 00167 END DO 00168 END IF 00169 WORK( 2*N+I ) = TMP 00170 END DO 00171 END IF 00172 * 00173 * Estimate the norm of inv(op(A)). 00174 * 00175 AINVNM = 0.0 00176 00177 KASE = 0 00178 10 CONTINUE 00179 CALL SLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE ) 00180 IF( KASE.NE.0 ) THEN 00181 IF( KASE.EQ.2 ) THEN 00182 * 00183 * Multiply by R. 00184 * 00185 DO I = 1, N 00186 WORK(I) = WORK(I) * WORK(2*N+I) 00187 END DO 00188 00189 IF (NOTRANS) THEN 00190 CALL SGETRS( 'No transpose', N, 1, AF, LDAF, IPIV, 00191 $ WORK, N, INFO ) 00192 ELSE 00193 CALL SGETRS( 'Transpose', N, 1, AF, LDAF, IPIV, 00194 $ WORK, N, INFO ) 00195 END IF 00196 * 00197 * Multiply by inv(C). 00198 * 00199 IF ( CMODE .EQ. 1 ) THEN 00200 DO I = 1, N 00201 WORK( I ) = WORK( I ) / C( I ) 00202 END DO 00203 ELSE IF ( CMODE .EQ. -1 ) THEN 00204 DO I = 1, N 00205 WORK( I ) = WORK( I ) * C( I ) 00206 END DO 00207 END IF 00208 ELSE 00209 * 00210 * Multiply by inv(C**T). 00211 * 00212 IF ( CMODE .EQ. 1 ) THEN 00213 DO I = 1, N 00214 WORK( I ) = WORK( I ) / C( I ) 00215 END DO 00216 ELSE IF ( CMODE .EQ. -1 ) THEN 00217 DO I = 1, N 00218 WORK( I ) = WORK( I ) * C( I ) 00219 END DO 00220 END IF 00221 00222 IF (NOTRANS) THEN 00223 CALL SGETRS( 'Transpose', N, 1, AF, LDAF, IPIV, 00224 $ WORK, N, INFO ) 00225 ELSE 00226 CALL SGETRS( 'No transpose', N, 1, AF, LDAF, IPIV, 00227 $ WORK, N, INFO ) 00228 END IF 00229 * 00230 * Multiply by R. 00231 * 00232 DO I = 1, N 00233 WORK( I ) = WORK( I ) * WORK( 2*N+I ) 00234 END DO 00235 END IF 00236 GO TO 10 00237 END IF 00238 * 00239 * Compute the estimate of the reciprocal condition number. 00240 * 00241 IF( AINVNM .NE. 0.0 ) 00242 $ SLA_GERCOND = ( 1.0 / AINVNM ) 00243 * 00244 RETURN 00245 * 00246 END
1.7.3 
