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