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