 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 DOUBLE PRECISION FUNCTION ZLA_SYRCOND_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_SYRCOND_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 ZSYTRF. 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 ZSYTRF. 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 00080 DOUBLE PRECISION AINVNM, ANORM, TMP 00081 INTEGER I, J 00082 LOGICAL UP 00083 COMPLEX*16 ZDUM 00084 * .. 00085 * .. Local Arrays .. 00086 INTEGER ISAVE( 3 ) 00087 * .. 00088 * .. External Functions .. 00089 LOGICAL LSAME 00090 EXTERNAL LSAME 00091 * .. 00092 * .. External Subroutines .. 00093 EXTERNAL ZLACN2, ZSYTRS, XERBLA 00094 * .. 00095 * .. Intrinsic Functions .. 00096 INTRINSIC ABS, MAX 00097 * .. 00098 * .. Statement Functions .. 00099 DOUBLE PRECISION CABS1 00100 * .. 00101 * .. Statement Function Definitions .. 00102 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 00103 * .. 00104 * .. Executable Statements .. 00105 * 00106 ZLA_SYRCOND_C = 0.0D+0 00107 * 00108 INFO = 0 00109 IF( N.LT.0 ) THEN 00110 INFO = -2 00111 END IF 00112 IF( INFO.NE.0 ) THEN 00113 CALL XERBLA( 'ZLA_SYRCOND_C', -INFO ) 00114 RETURN 00115 END IF 00116 UP = .FALSE. 00117 IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE. 00118 * 00119 * Compute norm of op(A)*op2(C). 00120 * 00121 ANORM = 0.0D+0 00122 IF ( UP ) THEN 00123 DO I = 1, N 00124 TMP = 0.0D+0 00125 IF ( CAPPLY ) THEN 00126 DO J = 1, I 00127 TMP = TMP + CABS1( A( J, I ) ) / C( J ) 00128 END DO 00129 DO J = I+1, N 00130 TMP = TMP + CABS1( A( I, J ) ) / C( J ) 00131 END DO 00132 ELSE 00133 DO J = 1, I 00134 TMP = TMP + CABS1( A( J, I ) ) 00135 END DO 00136 DO J = I+1, N 00137 TMP = TMP + CABS1( A( I, J ) ) 00138 END DO 00139 END IF 00140 RWORK( I ) = TMP 00141 ANORM = MAX( ANORM, TMP ) 00142 END DO 00143 ELSE 00144 DO I = 1, N 00145 TMP = 0.0D+0 00146 IF ( CAPPLY ) THEN 00147 DO J = 1, I 00148 TMP = TMP + CABS1( A( I, J ) ) / C( J ) 00149 END DO 00150 DO J = I+1, N 00151 TMP = TMP + CABS1( A( J, I ) ) / C( J ) 00152 END DO 00153 ELSE 00154 DO J = 1, I 00155 TMP = TMP + CABS1( A( I, J ) ) 00156 END DO 00157 DO J = I+1, N 00158 TMP = TMP + CABS1( A( J, I ) ) 00159 END DO 00160 END IF 00161 RWORK( I ) = TMP 00162 ANORM = MAX( ANORM, TMP ) 00163 END DO 00164 END IF 00165 * 00166 * Quick return if possible. 00167 * 00168 IF( N.EQ.0 ) THEN 00169 ZLA_SYRCOND_C = 1.0D+0 00170 RETURN 00171 ELSE IF( ANORM .EQ. 0.0D+0 ) THEN 00172 RETURN 00173 END IF 00174 * 00175 * Estimate the norm of inv(op(A)). 00176 * 00177 AINVNM = 0.0D+0 00178 * 00179 KASE = 0 00180 10 CONTINUE 00181 CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE ) 00182 IF( KASE.NE.0 ) THEN 00183 IF( KASE.EQ.2 ) THEN 00184 * 00185 * Multiply by R. 00186 * 00187 DO I = 1, N 00188 WORK( I ) = WORK( I ) * RWORK( I ) 00189 END DO 00190 * 00191 IF ( UP ) THEN 00192 CALL ZSYTRS( 'U', N, 1, AF, LDAF, IPIV, 00193 $ WORK, N, INFO ) 00194 ELSE 00195 CALL ZSYTRS( 'L', N, 1, AF, LDAF, IPIV, 00196 $ WORK, N, INFO ) 00197 ENDIF 00198 * 00199 * Multiply by inv(C). 00200 * 00201 IF ( CAPPLY ) THEN 00202 DO I = 1, N 00203 WORK( I ) = WORK( I ) * C( I ) 00204 END DO 00205 END IF 00206 ELSE 00207 * 00208 * Multiply by inv(C**T). 00209 * 00210 IF ( CAPPLY ) THEN 00211 DO I = 1, N 00212 WORK( I ) = WORK( I ) * C( I ) 00213 END DO 00214 END IF 00215 * 00216 IF ( UP ) THEN 00217 CALL ZSYTRS( 'U', N, 1, AF, LDAF, IPIV, 00218 $ WORK, N, INFO ) 00219 ELSE 00220 CALL ZSYTRS( 'L', N, 1, AF, LDAF, IPIV, 00221 $ WORK, N, INFO ) 00222 END IF 00223 * 00224 * Multiply by R. 00225 * 00226 DO I = 1, N 00227 WORK( I ) = WORK( I ) * RWORK( I ) 00228 END DO 00229 END IF 00230 GO TO 10 00231 END IF 00232 * 00233 * Compute the estimate of the reciprocal condition number. 00234 * 00235 IF( AINVNM .NE. 0.0D+0 ) 00236 $ ZLA_SYRCOND_C = 1.0D+0 / AINVNM 00237 * 00238 RETURN 00239 * 00240 END
1.7.3 
