 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 DOUBLE PRECISION FUNCTION ZLA_GERCOND_C( TRANS, N, A, LDA, AF, 00002 $ LDAF, IPIV, C, CAPPLY, 00003 $ 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 Aguments .. 00016 CHARACTER TRANS 00017 LOGICAL CAPPLY 00018 INTEGER N, LDA, LDAF, INFO 00019 * .. 00020 * .. Array Arguments .. 00021 INTEGER IPIV( * ) 00022 COMPLEX*16 A( LDA, * ), AF( LDAF, * ), WORK( * ) 00023 DOUBLE PRECISION C( * ), RWORK( * ) 00024 * .. 00025 * 00026 * Purpose 00027 * ======= 00028 * 00029 * ZLA_GERCOND_C computes the infinity norm condition number of 00030 * op(A) * inv(diag(C)) where C is a DOUBLE PRECISION 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 * A (input) COMPLEX*16 array, dimension (LDA,N) 00046 * On entry, the N-by-N matrix A 00047 * 00048 * LDA (input) INTEGER 00049 * The leading dimension of the array A. LDA >= max(1,N). 00050 * 00051 * AF (input) COMPLEX*16 array, dimension (LDAF,N) 00052 * The factors L and U from the factorization 00053 * A = P*L*U as computed by ZGETRF. 00054 * 00055 * LDAF (input) INTEGER 00056 * The leading dimension of the array AF. LDAF >= max(1,N). 00057 * 00058 * IPIV (input) INTEGER array, dimension (N) 00059 * The pivot indices from the factorization A = P*L*U 00060 * as computed by ZGETRF; row i of the matrix was interchanged 00061 * with row IPIV(i). 00062 * 00063 * C (input) DOUBLE PRECISION array, dimension (N) 00064 * The vector C in the formula op(A) * inv(diag(C)). 00065 * 00066 * CAPPLY (input) LOGICAL 00067 * If .TRUE. then access the vector C in the formula above. 00068 * 00069 * INFO (output) INTEGER 00070 * = 0: Successful exit. 00071 * i > 0: The ith argument is invalid. 00072 * 00073 * WORK (input) COMPLEX*16 array, dimension (2*N). 00074 * Workspace. 00075 * 00076 * RWORK (input) DOUBLE PRECISION array, dimension (N). 00077 * Workspace. 00078 * 00079 * ===================================================================== 00080 * 00081 * .. Local Scalars .. 00082 LOGICAL NOTRANS 00083 INTEGER KASE, I, J 00084 DOUBLE PRECISION AINVNM, ANORM, TMP 00085 COMPLEX*16 ZDUM 00086 * .. 00087 * .. Local Arrays .. 00088 INTEGER ISAVE( 3 ) 00089 * .. 00090 * .. External Functions .. 00091 LOGICAL LSAME 00092 EXTERNAL LSAME 00093 * .. 00094 * .. External Subroutines .. 00095 EXTERNAL ZLACN2, ZGETRS, XERBLA 00096 * .. 00097 * .. Intrinsic Functions .. 00098 INTRINSIC ABS, MAX, REAL, DIMAG 00099 * .. 00100 * .. Statement Functions .. 00101 DOUBLE PRECISION CABS1 00102 * .. 00103 * .. Statement Function Definitions .. 00104 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 00105 * .. 00106 * .. Executable Statements .. 00107 ZLA_GERCOND_C = 0.0D+0 00108 * 00109 INFO = 0 00110 NOTRANS = LSAME( TRANS, 'N' ) 00111 IF ( .NOT. NOTRANS .AND. .NOT. LSAME( TRANS, 'T' ) .AND. .NOT. 00112 $ LSAME( TRANS, 'C' ) ) THEN 00113 ELSE IF( N.LT.0 ) THEN 00114 INFO = -2 00115 END IF 00116 IF( INFO.NE.0 ) THEN 00117 CALL XERBLA( 'ZLA_GERCOND_C', -INFO ) 00118 RETURN 00119 END IF 00120 * 00121 * Compute norm of op(A)*op2(C). 00122 * 00123 ANORM = 0.0D+0 00124 IF ( NOTRANS ) THEN 00125 DO I = 1, N 00126 TMP = 0.0D+0 00127 IF ( CAPPLY ) THEN 00128 DO J = 1, N 00129 TMP = TMP + CABS1( A( I, J ) ) / C( J ) 00130 END DO 00131 ELSE 00132 DO J = 1, N 00133 TMP = TMP + CABS1( A( I, J ) ) 00134 END DO 00135 END IF 00136 RWORK( I ) = TMP 00137 ANORM = MAX( ANORM, TMP ) 00138 END DO 00139 ELSE 00140 DO I = 1, N 00141 TMP = 0.0D+0 00142 IF ( CAPPLY ) THEN 00143 DO J = 1, N 00144 TMP = TMP + CABS1( A( J, I ) ) / C( J ) 00145 END DO 00146 ELSE 00147 DO J = 1, N 00148 TMP = TMP + CABS1( A( J, I ) ) 00149 END DO 00150 END IF 00151 RWORK( I ) = TMP 00152 ANORM = MAX( ANORM, TMP ) 00153 END DO 00154 END IF 00155 * 00156 * Quick return if possible. 00157 * 00158 IF( N.EQ.0 ) THEN 00159 ZLA_GERCOND_C = 1.0D+0 00160 RETURN 00161 ELSE IF( ANORM .EQ. 0.0D+0 ) THEN 00162 RETURN 00163 END IF 00164 * 00165 * Estimate the norm of inv(op(A)). 00166 * 00167 AINVNM = 0.0D+0 00168 * 00169 KASE = 0 00170 10 CONTINUE 00171 CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE ) 00172 IF( KASE.NE.0 ) THEN 00173 IF( KASE.EQ.2 ) THEN 00174 * 00175 * Multiply by R. 00176 * 00177 DO I = 1, N 00178 WORK( I ) = WORK( I ) * RWORK( I ) 00179 END DO 00180 * 00181 IF (NOTRANS) THEN 00182 CALL ZGETRS( 'No transpose', N, 1, AF, LDAF, IPIV, 00183 $ WORK, N, INFO ) 00184 ELSE 00185 CALL ZGETRS( 'Conjugate transpose', N, 1, AF, LDAF, IPIV, 00186 $ WORK, N, INFO ) 00187 ENDIF 00188 * 00189 * Multiply by inv(C). 00190 * 00191 IF ( CAPPLY ) THEN 00192 DO I = 1, N 00193 WORK( I ) = WORK( I ) * C( I ) 00194 END DO 00195 END IF 00196 ELSE 00197 * 00198 * Multiply by inv(C**H). 00199 * 00200 IF ( CAPPLY ) THEN 00201 DO I = 1, N 00202 WORK( I ) = WORK( I ) * C( I ) 00203 END DO 00204 END IF 00205 * 00206 IF ( NOTRANS ) THEN 00207 CALL ZGETRS( 'Conjugate transpose', N, 1, AF, LDAF, IPIV, 00208 $ WORK, N, INFO ) 00209 ELSE 00210 CALL ZGETRS( 'No transpose', N, 1, AF, LDAF, IPIV, 00211 $ WORK, N, INFO ) 00212 END IF 00213 * 00214 * Multiply by R. 00215 * 00216 DO I = 1, N 00217 WORK( I ) = WORK( I ) * RWORK( I ) 00218 END DO 00219 END IF 00220 GO TO 10 00221 END IF 00222 * 00223 * Compute the estimate of the reciprocal condition number. 00224 * 00225 IF( AINVNM .NE. 0.0D+0 ) 00226 $ ZLA_GERCOND_C = 1.0D+0 / AINVNM 00227 * 00228 RETURN 00229 * 00230 END
1.7.3 
