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