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