|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 DOUBLE PRECISION FUNCTION DLA_GBRCOND( TRANS, N, KL, KU, AB, LDAB, 00002 $ AFB, LDAFB, IPIV, CMODE, C, 00003 $ INFO, WORK, IWORK ) 00004 * 00005 * -- LAPACK routine (version 3.2.2) -- 00006 * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and -- 00007 * -- Jason Riedy of Univ. of California Berkeley. -- 00008 * -- June 2010 -- 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, LDAB, LDAFB, INFO, KL, KU, CMODE 00018 * .. 00019 * .. Array Arguments .. 00020 INTEGER IWORK( * ), IPIV( * ) 00021 DOUBLE PRECISION AB( LDAB, * ), AFB( LDAFB, * ), WORK( * ), 00022 $ C( * ) 00023 * .. 00024 * 00025 * Purpose 00026 * ======= 00027 * 00028 * DLA_GBRCOND 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 * KL (input) INTEGER 00052 * The number of subdiagonals within the band of A. KL >= 0. 00053 * 00054 * KU (input) INTEGER 00055 * The number of superdiagonals within the band of A. KU >= 0. 00056 * 00057 * AB (input) DOUBLE PRECISION array, dimension (LDAB,N) 00058 * On entry, the matrix A in band storage, in rows 1 to KL+KU+1. 00059 * The j-th column of A is stored in the j-th column of the 00060 * array AB as follows: 00061 * AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl) 00062 * 00063 * LDAB (input) INTEGER 00064 * The leading dimension of the array AB. LDAB >= KL+KU+1. 00065 * 00066 * AFB (input) DOUBLE PRECISION array, dimension (LDAFB,N) 00067 * Details of the LU factorization of the band matrix A, as 00068 * computed by DGBTRF. U is stored as an upper triangular 00069 * band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, 00070 * and the multipliers used during the factorization are stored 00071 * in rows KL+KU+2 to 2*KL+KU+1. 00072 * 00073 * LDAFB (input) INTEGER 00074 * The leading dimension of the array AFB. LDAFB >= 2*KL+KU+1. 00075 * 00076 * IPIV (input) INTEGER array, dimension (N) 00077 * The pivot indices from the factorization A = P*L*U 00078 * as computed by DGBTRF; row i of the matrix was interchanged 00079 * with row IPIV(i). 00080 * 00081 * CMODE (input) INTEGER 00082 * Determines op2(C) in the formula op(A) * op2(C) as follows: 00083 * CMODE = 1 op2(C) = C 00084 * CMODE = 0 op2(C) = I 00085 * CMODE = -1 op2(C) = inv(C) 00086 * 00087 * C (input) DOUBLE PRECISION array, dimension (N) 00088 * The vector C in the formula op(A) * op2(C). 00089 * 00090 * INFO (output) INTEGER 00091 * = 0: Successful exit. 00092 * i > 0: The ith argument is invalid. 00093 * 00094 * WORK (input) DOUBLE PRECISION array, dimension (5*N). 00095 * Workspace. 00096 * 00097 * IWORK (input) INTEGER array, dimension (N). 00098 * Workspace. 00099 * 00100 * ===================================================================== 00101 * 00102 * .. Local Scalars .. 00103 LOGICAL NOTRANS 00104 INTEGER KASE, I, J, KD, KE 00105 DOUBLE PRECISION AINVNM, TMP 00106 * .. 00107 * .. Local Arrays .. 00108 INTEGER ISAVE( 3 ) 00109 * .. 00110 * .. External Functions .. 00111 LOGICAL LSAME 00112 EXTERNAL LSAME 00113 * .. 00114 * .. External Subroutines .. 00115 EXTERNAL DLACN2, DGBTRS, XERBLA 00116 * .. 00117 * .. Intrinsic Functions .. 00118 INTRINSIC ABS, MAX 00119 * .. 00120 * .. Executable Statements .. 00121 * 00122 DLA_GBRCOND = 0.0D+0 00123 * 00124 INFO = 0 00125 NOTRANS = LSAME( TRANS, 'N' ) 00126 IF ( .NOT. NOTRANS .AND. .NOT. LSAME(TRANS, 'T') 00127 $ .AND. .NOT. LSAME(TRANS, 'C') ) THEN 00128 INFO = -1 00129 ELSE IF( N.LT.0 ) THEN 00130 INFO = -2 00131 ELSE IF( KL.LT.0 .OR. KL.GT.N-1 ) THEN 00132 INFO = -3 00133 ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN 00134 INFO = -4 00135 ELSE IF( LDAB.LT.KL+KU+1 ) THEN 00136 INFO = -6 00137 ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN 00138 INFO = -8 00139 END IF 00140 IF( INFO.NE.0 ) THEN 00141 CALL XERBLA( 'DLA_GBRCOND', -INFO ) 00142 RETURN 00143 END IF 00144 IF( N.EQ.0 ) THEN 00145 DLA_GBRCOND = 1.0D+0 00146 RETURN 00147 END IF 00148 * 00149 * Compute the equilibration matrix R such that 00150 * inv(R)*A*C has unit 1-norm. 00151 * 00152 KD = KU + 1 00153 KE = KL + 1 00154 IF ( NOTRANS ) THEN 00155 DO I = 1, N 00156 TMP = 0.0D+0 00157 IF ( CMODE .EQ. 1 ) THEN 00158 DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 00159 TMP = TMP + ABS( AB( KD+I-J, J ) * C( J ) ) 00160 END DO 00161 ELSE IF ( CMODE .EQ. 0 ) THEN 00162 DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 00163 TMP = TMP + ABS( AB( KD+I-J, J ) ) 00164 END DO 00165 ELSE 00166 DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 00167 TMP = TMP + ABS( AB( KD+I-J, J ) / C( J ) ) 00168 END DO 00169 END IF 00170 WORK( 2*N+I ) = TMP 00171 END DO 00172 ELSE 00173 DO I = 1, N 00174 TMP = 0.0D+0 00175 IF ( CMODE .EQ. 1 ) THEN 00176 DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 00177 TMP = TMP + ABS( AB( KE-I+J, I ) * C( J ) ) 00178 END DO 00179 ELSE IF ( CMODE .EQ. 0 ) THEN 00180 DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 00181 TMP = TMP + ABS( AB( KE-I+J, I ) ) 00182 END DO 00183 ELSE 00184 DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 00185 TMP = TMP + ABS( AB( KE-I+J, I ) / C( J ) ) 00186 END DO 00187 END IF 00188 WORK( 2*N+I ) = TMP 00189 END DO 00190 END IF 00191 * 00192 * Estimate the norm of inv(op(A)). 00193 * 00194 AINVNM = 0.0D+0 00195 00196 KASE = 0 00197 10 CONTINUE 00198 CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE ) 00199 IF( KASE.NE.0 ) THEN 00200 IF( KASE.EQ.2 ) THEN 00201 * 00202 * Multiply by R. 00203 * 00204 DO I = 1, N 00205 WORK( I ) = WORK( I ) * WORK( 2*N+I ) 00206 END DO 00207 00208 IF ( NOTRANS ) THEN 00209 CALL DGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB, 00210 $ IPIV, WORK, N, INFO ) 00211 ELSE 00212 CALL DGBTRS( 'Transpose', N, KL, KU, 1, AFB, LDAFB, IPIV, 00213 $ WORK, N, INFO ) 00214 END IF 00215 * 00216 * Multiply by inv(C). 00217 * 00218 IF ( CMODE .EQ. 1 ) THEN 00219 DO I = 1, N 00220 WORK( I ) = WORK( I ) / C( I ) 00221 END DO 00222 ELSE IF ( CMODE .EQ. -1 ) THEN 00223 DO I = 1, N 00224 WORK( I ) = WORK( I ) * C( I ) 00225 END DO 00226 END IF 00227 ELSE 00228 * 00229 * Multiply by inv(C'). 00230 * 00231 IF ( CMODE .EQ. 1 ) THEN 00232 DO I = 1, N 00233 WORK( I ) = WORK( I ) / C( I ) 00234 END DO 00235 ELSE IF ( CMODE .EQ. -1 ) THEN 00236 DO I = 1, N 00237 WORK( I ) = WORK( I ) * C( I ) 00238 END DO 00239 END IF 00240 00241 IF ( NOTRANS ) THEN 00242 CALL DGBTRS( 'Transpose', N, KL, KU, 1, AFB, LDAFB, IPIV, 00243 $ WORK, N, INFO ) 00244 ELSE 00245 CALL DGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB, 00246 $ IPIV, WORK, N, INFO ) 00247 END IF 00248 * 00249 * Multiply by R. 00250 * 00251 DO I = 1, N 00252 WORK( I ) = WORK( I ) * WORK( 2*N+I ) 00253 END DO 00254 END IF 00255 GO TO 10 00256 END IF 00257 * 00258 * Compute the estimate of the reciprocal condition number. 00259 * 00260 IF( AINVNM .NE. 0.0D+0 ) 00261 $ DLA_GBRCOND = ( 1.0D+0 / AINVNM ) 00262 * 00263 RETURN 00264 * 00265 END
1.7.3 
