 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 DOUBLE PRECISION FUNCTION ZLA_HERCOND_C( UPLO, 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 Arguments .. 00016 CHARACTER UPLO 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_HERCOND_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 * UPLO (input) CHARACTER*1 00036 * = 'U': Upper triangle of A is stored; 00037 * = 'L': Lower triangle of A is stored. 00038 * 00039 * N (input) INTEGER 00040 * The number of linear equations, i.e., the order of the 00041 * matrix A. N >= 0. 00042 * 00043 * A (input) COMPLEX*16 array, dimension (LDA,N) 00044 * On entry, the N-by-N matrix A 00045 * 00046 * LDA (input) INTEGER 00047 * The leading dimension of the array A. LDA >= max(1,N). 00048 * 00049 * AF (input) COMPLEX*16 array, dimension (LDAF,N) 00050 * The block diagonal matrix D and the multipliers used to 00051 * obtain the factor U or L as computed by ZHETRF. 00052 * 00053 * LDAF (input) INTEGER 00054 * The leading dimension of the array AF. LDAF >= max(1,N). 00055 * 00056 * IPIV (input) INTEGER array, dimension (N) 00057 * Details of the interchanges and the block structure of D 00058 * as determined by CHETRF. 00059 * 00060 * C (input) DOUBLE PRECISION array, dimension (N) 00061 * The vector C in the formula op(A) * inv(diag(C)). 00062 * 00063 * CAPPLY (input) LOGICAL 00064 * If .TRUE. then access the vector C in the formula above. 00065 * 00066 * INFO (output) INTEGER 00067 * = 0: Successful exit. 00068 * i > 0: The ith argument is invalid. 00069 * 00070 * WORK (input) COMPLEX*16 array, dimension (2*N). 00071 * Workspace. 00072 * 00073 * RWORK (input) DOUBLE PRECISION array, dimension (N). 00074 * Workspace. 00075 * 00076 * ===================================================================== 00077 * 00078 * .. Local Scalars .. 00079 INTEGER KASE, I, J 00080 DOUBLE PRECISION AINVNM, ANORM, TMP 00081 LOGICAL UP 00082 COMPLEX*16 ZDUM 00083 * .. 00084 * .. Local Arrays .. 00085 INTEGER ISAVE( 3 ) 00086 * .. 00087 * .. External Functions .. 00088 LOGICAL LSAME 00089 EXTERNAL LSAME 00090 * .. 00091 * .. External Subroutines .. 00092 EXTERNAL ZLACN2, ZHETRS, XERBLA 00093 * .. 00094 * .. Intrinsic Functions .. 00095 INTRINSIC ABS, MAX 00096 * .. 00097 * .. Statement Functions .. 00098 DOUBLE PRECISION CABS1 00099 * .. 00100 * .. Statement Function Definitions .. 00101 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 00102 * .. 00103 * .. Executable Statements .. 00104 * 00105 ZLA_HERCOND_C = 0.0D+0 00106 * 00107 INFO = 0 00108 IF( N.LT.0 ) THEN 00109 INFO = -2 00110 END IF 00111 IF( INFO.NE.0 ) THEN 00112 CALL XERBLA( 'ZLA_HERCOND_C', -INFO ) 00113 RETURN 00114 END IF 00115 UP = .FALSE. 00116 IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE. 00117 * 00118 * Compute norm of op(A)*op2(C). 00119 * 00120 ANORM = 0.0D+0 00121 IF ( UP ) THEN 00122 DO I = 1, N 00123 TMP = 0.0D+0 00124 IF ( CAPPLY ) THEN 00125 DO J = 1, I 00126 TMP = TMP + CABS1( A( J, I ) ) / C( J ) 00127 END DO 00128 DO J = I+1, N 00129 TMP = TMP + CABS1( A( I, J ) ) / C( J ) 00130 END DO 00131 ELSE 00132 DO J = 1, I 00133 TMP = TMP + CABS1( A( J, I ) ) 00134 END DO 00135 DO J = I+1, N 00136 TMP = TMP + CABS1( A( I, J ) ) 00137 END DO 00138 END IF 00139 RWORK( I ) = TMP 00140 ANORM = MAX( ANORM, TMP ) 00141 END DO 00142 ELSE 00143 DO I = 1, N 00144 TMP = 0.0D+0 00145 IF ( CAPPLY ) THEN 00146 DO J = 1, I 00147 TMP = TMP + CABS1( A( I, J ) ) / C( J ) 00148 END DO 00149 DO J = I+1, N 00150 TMP = TMP + CABS1( A( J, I ) ) / C( J ) 00151 END DO 00152 ELSE 00153 DO J = 1, I 00154 TMP = TMP + CABS1( A( I, J ) ) 00155 END DO 00156 DO J = I+1, N 00157 TMP = TMP + CABS1( A( J, I ) ) 00158 END DO 00159 END IF 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_HERCOND_C = 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 ( UP ) THEN 00191 CALL ZHETRS( 'U', N, 1, AF, LDAF, IPIV, 00192 $ WORK, N, INFO ) 00193 ELSE 00194 CALL ZHETRS( 'L', N, 1, AF, LDAF, IPIV, 00195 $ WORK, N, INFO ) 00196 ENDIF 00197 * 00198 * Multiply by inv(C). 00199 * 00200 IF ( CAPPLY ) THEN 00201 DO I = 1, N 00202 WORK( I ) = WORK( I ) * C( I ) 00203 END DO 00204 END IF 00205 ELSE 00206 * 00207 * Multiply by inv(C**H). 00208 * 00209 IF ( CAPPLY ) THEN 00210 DO I = 1, N 00211 WORK( I ) = WORK( I ) * C( I ) 00212 END DO 00213 END IF 00214 * 00215 IF ( UP ) THEN 00216 CALL ZHETRS( 'U', N, 1, AF, LDAF, IPIV, 00217 $ WORK, N, INFO ) 00218 ELSE 00219 CALL ZHETRS( 'L', N, 1, AF, LDAF, IPIV, 00220 $ WORK, N, INFO ) 00221 END IF 00222 * 00223 * Multiply by R. 00224 * 00225 DO I = 1, N 00226 WORK( I ) = WORK( I ) * RWORK( I ) 00227 END DO 00228 END IF 00229 GO TO 10 00230 END IF 00231 * 00232 * Compute the estimate of the reciprocal condition number. 00233 * 00234 IF( AINVNM .NE. 0.0D+0 ) 00235 $ ZLA_HERCOND_C = 1.0D+0 / AINVNM 00236 * 00237 RETURN 00238 * 00239 END
1.7.3 
