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