 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 REAL FUNCTION CLA_GBRCOND_C( TRANS, N, KL, KU, AB, LDAB, AFB, 00002 $ LDAFB, IPIV, C, CAPPLY, INFO, WORK, 00003 $ RWORK ) 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 LOGICAL CAPPLY 00018 INTEGER N, KL, KU, KD, KE, LDAB, LDAFB, INFO 00019 * .. 00020 * .. Array Arguments .. 00021 INTEGER IPIV( * ) 00022 COMPLEX AB( LDAB, * ), AFB( LDAFB, * ), WORK( * ) 00023 REAL C( * ), RWORK( * ) 00024 * .. 00025 * 00026 * Purpose 00027 * ======= 00028 * 00029 * CLA_GBRCOND_C Computes the infinity norm condition number of 00030 * op(A) * inv(diag(C)) where C is a REAL vector. 00031 * 00032 * Arguments 00033 * ========= 00034 * 00035 * TRANS (input) CHARACTER*1 00036 * Specifies the form of the system of equations: 00037 * = 'N': A * X = B (No transpose) 00038 * = 'T': A**T * X = B (Transpose) 00039 * = 'C': A**H * X = B (Conjugate Transpose = Transpose) 00040 * 00041 * N (input) INTEGER 00042 * The number of linear equations, i.e., the order of the 00043 * matrix A. N >= 0. 00044 * 00045 * KL (input) INTEGER 00046 * The number of subdiagonals within the band of A. KL >= 0. 00047 * 00048 * KU (input) INTEGER 00049 * The number of superdiagonals within the band of A. KU >= 0. 00050 * 00051 * AB (input) COMPLEX array, dimension (LDAB,N) 00052 * On entry, the matrix A in band storage, in rows 1 to KL+KU+1. 00053 * The j-th column of A is stored in the j-th column of the 00054 * array AB as follows: 00055 * AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl) 00056 * 00057 * LDAB (input) INTEGER 00058 * The leading dimension of the array AB. LDAB >= KL+KU+1. 00059 * 00060 * AFB (input) COMPLEX array, dimension (LDAFB,N) 00061 * Details of the LU factorization of the band matrix A, as 00062 * computed by CGBTRF. U is stored as an upper triangular 00063 * band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, 00064 * and the multipliers used during the factorization are stored 00065 * in rows KL+KU+2 to 2*KL+KU+1. 00066 * 00067 * LDAFB (input) INTEGER 00068 * The leading dimension of the array AFB. LDAFB >= 2*KL+KU+1. 00069 * 00070 * IPIV (input) INTEGER array, dimension (N) 00071 * The pivot indices from the factorization A = P*L*U 00072 * as computed by CGBTRF; row i of the matrix was interchanged 00073 * with row IPIV(i). 00074 * 00075 * C (input) REAL array, dimension (N) 00076 * The vector C in the formula op(A) * inv(diag(C)). 00077 * 00078 * CAPPLY (input) LOGICAL 00079 * If .TRUE. then access the vector C in the formula above. 00080 * 00081 * INFO (output) INTEGER 00082 * = 0: Successful exit. 00083 * i > 0: The ith argument is invalid. 00084 * 00085 * WORK (input) COMPLEX array, dimension (2*N). 00086 * Workspace. 00087 * 00088 * RWORK (input) REAL array, dimension (N). 00089 * Workspace. 00090 * 00091 * ===================================================================== 00092 * 00093 * .. Local Scalars .. 00094 LOGICAL NOTRANS 00095 INTEGER KASE, I, J 00096 REAL AINVNM, ANORM, TMP 00097 COMPLEX ZDUM 00098 * .. 00099 * .. Local Arrays .. 00100 INTEGER ISAVE( 3 ) 00101 * .. 00102 * .. External Functions .. 00103 LOGICAL LSAME 00104 EXTERNAL LSAME 00105 * .. 00106 * .. External Subroutines .. 00107 EXTERNAL CLACN2, CGBTRS, XERBLA 00108 * .. 00109 * .. Intrinsic Functions .. 00110 INTRINSIC ABS, MAX 00111 * .. 00112 * .. Statement Functions .. 00113 REAL CABS1 00114 * .. 00115 * .. Statement Function Definitions .. 00116 CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) ) 00117 * .. 00118 * .. Executable Statements .. 00119 CLA_GBRCOND_C = 0.0E+0 00120 * 00121 INFO = 0 00122 NOTRANS = LSAME( TRANS, 'N' ) 00123 IF ( .NOT. NOTRANS .AND. .NOT. LSAME( TRANS, 'T' ) .AND. .NOT. 00124 $ LSAME( TRANS, 'C' ) ) THEN 00125 INFO = -1 00126 ELSE IF( N.LT.0 ) THEN 00127 INFO = -2 00128 ELSE IF( KL.LT.0 .OR. KL.GT.N-1 ) THEN 00129 INFO = -3 00130 ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN 00131 INFO = -4 00132 ELSE IF( LDAB.LT.KL+KU+1 ) THEN 00133 INFO = -6 00134 ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN 00135 INFO = -8 00136 END IF 00137 IF( INFO.NE.0 ) THEN 00138 CALL XERBLA( 'CLA_GBRCOND_C', -INFO ) 00139 RETURN 00140 END IF 00141 * 00142 * Compute norm of op(A)*op2(C). 00143 * 00144 ANORM = 0.0E+0 00145 KD = KU + 1 00146 KE = KL + 1 00147 IF ( NOTRANS ) THEN 00148 DO I = 1, N 00149 TMP = 0.0E+0 00150 IF ( CAPPLY ) THEN 00151 DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 00152 TMP = TMP + CABS1( AB( KD+I-J, J ) ) / C( J ) 00153 END DO 00154 ELSE 00155 DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 00156 TMP = TMP + CABS1( AB( KD+I-J, J ) ) 00157 END DO 00158 END IF 00159 RWORK( I ) = TMP 00160 ANORM = MAX( ANORM, TMP ) 00161 END DO 00162 ELSE 00163 DO I = 1, N 00164 TMP = 0.0E+0 00165 IF ( CAPPLY ) THEN 00166 DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 00167 TMP = TMP + CABS1( AB( KE-I+J, I ) ) / C( J ) 00168 END DO 00169 ELSE 00170 DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 00171 TMP = TMP + CABS1( AB( KE-I+J, I ) ) 00172 END DO 00173 END IF 00174 RWORK( I ) = TMP 00175 ANORM = MAX( ANORM, TMP ) 00176 END DO 00177 END IF 00178 * 00179 * Quick return if possible. 00180 * 00181 IF( N.EQ.0 ) THEN 00182 CLA_GBRCOND_C = 1.0E+0 00183 RETURN 00184 ELSE IF( ANORM .EQ. 0.0E+0 ) THEN 00185 RETURN 00186 END IF 00187 * 00188 * Estimate the norm of inv(op(A)). 00189 * 00190 AINVNM = 0.0E+0 00191 * 00192 KASE = 0 00193 10 CONTINUE 00194 CALL CLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE ) 00195 IF( KASE.NE.0 ) THEN 00196 IF( KASE.EQ.2 ) THEN 00197 * 00198 * Multiply by R. 00199 * 00200 DO I = 1, N 00201 WORK( I ) = WORK( I ) * RWORK( I ) 00202 END DO 00203 * 00204 IF ( NOTRANS ) THEN 00205 CALL CGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB, 00206 $ IPIV, WORK, N, INFO ) 00207 ELSE 00208 CALL CGBTRS( 'Conjugate transpose', N, KL, KU, 1, AFB, 00209 $ LDAFB, IPIV, WORK, N, INFO ) 00210 ENDIF 00211 * 00212 * Multiply by inv(C). 00213 * 00214 IF ( CAPPLY ) THEN 00215 DO I = 1, N 00216 WORK( I ) = WORK( I ) * C( I ) 00217 END DO 00218 END IF 00219 ELSE 00220 * 00221 * Multiply by inv(C**H). 00222 * 00223 IF ( CAPPLY ) THEN 00224 DO I = 1, N 00225 WORK( I ) = WORK( I ) * C( I ) 00226 END DO 00227 END IF 00228 * 00229 IF ( NOTRANS ) THEN 00230 CALL CGBTRS( 'Conjugate transpose', N, KL, KU, 1, AFB, 00231 $ LDAFB, IPIV, WORK, N, INFO ) 00232 ELSE 00233 CALL CGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB, 00234 $ IPIV, WORK, N, INFO ) 00235 END IF 00236 * 00237 * Multiply by R. 00238 * 00239 DO I = 1, N 00240 WORK( I ) = WORK( I ) * RWORK( I ) 00241 END DO 00242 END IF 00243 GO TO 10 00244 END IF 00245 * 00246 * Compute the estimate of the reciprocal condition number. 00247 * 00248 IF( AINVNM .NE. 0.0E+0 ) 00249 $ CLA_GBRCOND_C = 1.0E+0 / AINVNM 00250 * 00251 RETURN 00252 * 00253 END
1.7.3 
