|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 DOUBLE PRECISION FUNCTION ZLA_GBRCOND_X( TRANS, N, KL, KU, AB, 00002 $ LDAB, AFB, LDAFB, IPIV, 00003 $ X, INFO, WORK, 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 INTEGER N, KL, KU, KD, KE, LDAB, LDAFB, INFO 00018 * .. 00019 * .. Array Arguments .. 00020 INTEGER IPIV( * ) 00021 COMPLEX*16 AB( LDAB, * ), AFB( LDAFB, * ), WORK( * ), 00022 $ X( * ) 00023 DOUBLE PRECISION RWORK( * ) 00024 * 00025 * 00026 * Purpose 00027 * ======= 00028 * 00029 * ZLA_GBRCOND_X Computes the infinity norm condition number of 00030 * op(A) * diag(X) where X is a COMPLEX*16 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*16 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*16 array, dimension (LDAFB,N) 00061 * Details of the LU factorization of the band matrix A, as 00062 * computed by ZGBTRF. 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 ZGBTRF; row i of the matrix was interchanged 00073 * with row IPIV(i). 00074 * 00075 * X (input) COMPLEX*16 array, dimension (N) 00076 * The vector X in the formula op(A) * diag(X). 00077 * 00078 * INFO (output) INTEGER 00079 * = 0: Successful exit. 00080 * i > 0: The ith argument is invalid. 00081 * 00082 * WORK (input) COMPLEX*16 array, dimension (2*N). 00083 * Workspace. 00084 * 00085 * RWORK (input) DOUBLE PRECISION array, dimension (N). 00086 * Workspace. 00087 * 00088 * ===================================================================== 00089 * 00090 * .. Local Scalars .. 00091 LOGICAL NOTRANS 00092 INTEGER KASE, I, J 00093 DOUBLE PRECISION AINVNM, ANORM, TMP 00094 COMPLEX*16 ZDUM 00095 * .. 00096 * .. Local Arrays .. 00097 INTEGER ISAVE( 3 ) 00098 * .. 00099 * .. External Functions .. 00100 LOGICAL LSAME 00101 EXTERNAL LSAME 00102 * .. 00103 * .. External Subroutines .. 00104 EXTERNAL ZLACN2, ZGBTRS, XERBLA 00105 * .. 00106 * .. Intrinsic Functions .. 00107 INTRINSIC ABS, MAX 00108 * .. 00109 * .. Statement Functions .. 00110 DOUBLE PRECISION CABS1 00111 * .. 00112 * .. Statement Function Definitions .. 00113 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 00114 * .. 00115 * .. Executable Statements .. 00116 * 00117 ZLA_GBRCOND_X = 0.0D+0 00118 * 00119 INFO = 0 00120 NOTRANS = LSAME( TRANS, 'N' ) 00121 IF ( .NOT. NOTRANS .AND. .NOT. LSAME(TRANS, 'T') .AND. .NOT. 00122 $ LSAME( TRANS, 'C' ) ) THEN 00123 INFO = -1 00124 ELSE IF( N.LT.0 ) THEN 00125 INFO = -2 00126 ELSE IF( KL.LT.0 .OR. KL.GT.N-1 ) THEN 00127 INFO = -3 00128 ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN 00129 INFO = -4 00130 ELSE IF( LDAB.LT.KL+KU+1 ) THEN 00131 INFO = -6 00132 ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN 00133 INFO = -8 00134 END IF 00135 IF( INFO.NE.0 ) THEN 00136 CALL XERBLA( 'ZLA_GBRCOND_X', -INFO ) 00137 RETURN 00138 END IF 00139 * 00140 * Compute norm of op(A)*op2(C). 00141 * 00142 KD = KU + 1 00143 KE = KL + 1 00144 ANORM = 0.0D+0 00145 IF ( NOTRANS ) THEN 00146 DO I = 1, N 00147 TMP = 0.0D+0 00148 DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 00149 TMP = TMP + CABS1( AB( KD+I-J, J) * X( J ) ) 00150 END DO 00151 RWORK( I ) = TMP 00152 ANORM = MAX( ANORM, TMP ) 00153 END DO 00154 ELSE 00155 DO I = 1, N 00156 TMP = 0.0D+0 00157 DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 00158 TMP = TMP + CABS1( AB( KE-I+J, I ) * X( J ) ) 00159 END DO 00160 RWORK( I ) = TMP 00161 ANORM = MAX( ANORM, TMP ) 00162 END DO 00163 END IF 00164 * 00165 * Quick return if possible. 00166 * 00167 IF( N.EQ.0 ) THEN 00168 ZLA_GBRCOND_X = 1.0D+0 00169 RETURN 00170 ELSE IF( ANORM .EQ. 0.0D+0 ) THEN 00171 RETURN 00172 END IF 00173 * 00174 * Estimate the norm of inv(op(A)). 00175 * 00176 AINVNM = 0.0D+0 00177 * 00178 KASE = 0 00179 10 CONTINUE 00180 CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE ) 00181 IF( KASE.NE.0 ) THEN 00182 IF( KASE.EQ.2 ) THEN 00183 * 00184 * Multiply by R. 00185 * 00186 DO I = 1, N 00187 WORK( I ) = WORK( I ) * RWORK( I ) 00188 END DO 00189 * 00190 IF ( NOTRANS ) THEN 00191 CALL ZGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB, 00192 $ IPIV, WORK, N, INFO ) 00193 ELSE 00194 CALL ZGBTRS( 'Conjugate transpose', N, KL, KU, 1, AFB, 00195 $ LDAFB, IPIV, WORK, N, INFO ) 00196 ENDIF 00197 * 00198 * Multiply by inv(X). 00199 * 00200 DO I = 1, N 00201 WORK( I ) = WORK( I ) / X( I ) 00202 END DO 00203 ELSE 00204 * 00205 * Multiply by inv(X'). 00206 * 00207 DO I = 1, N 00208 WORK( I ) = WORK( I ) / X( I ) 00209 END DO 00210 * 00211 IF ( NOTRANS ) THEN 00212 CALL ZGBTRS( 'Conjugate transpose', N, KL, KU, 1, AFB, 00213 $ LDAFB, IPIV, WORK, N, INFO ) 00214 ELSE 00215 CALL ZGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB, 00216 $ IPIV, WORK, N, INFO ) 00217 END IF 00218 * 00219 * Multiply by R. 00220 * 00221 DO I = 1, N 00222 WORK( I ) = WORK( I ) * RWORK( I ) 00223 END DO 00224 END IF 00225 GO TO 10 00226 END IF 00227 * 00228 * Compute the estimate of the reciprocal condition number. 00229 * 00230 IF( AINVNM .NE. 0.0D+0 ) 00231 $ ZLA_GBRCOND_X = 1.0D+0 / AINVNM 00232 * 00233 RETURN 00234 * 00235 END
1.7.3 
