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