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