![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLA_PORCOND 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DLA_PORCOND + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dla_porcond.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dla_porcond.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dla_porcond.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * DOUBLE PRECISION FUNCTION DLA_PORCOND( UPLO, N, A, LDA, AF, LDAF, 00022 * CMODE, C, INFO, WORK, 00023 * IWORK ) 00024 * 00025 * .. Scalar Arguments .. 00026 * CHARACTER UPLO 00027 * INTEGER N, LDA, LDAF, INFO, CMODE 00028 * DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), WORK( * ), 00029 * $ C( * ) 00030 * .. 00031 * .. Array Arguments .. 00032 * INTEGER IWORK( * ) 00033 * .. 00034 * 00035 * 00036 *> \par Purpose: 00037 * ============= 00038 *> 00039 *> \verbatim 00040 *> 00041 *> DLA_PORCOND Estimates the Skeel condition number of op(A) * op2(C) 00042 *> where op2 is determined by CMODE as follows 00043 *> CMODE = 1 op2(C) = C 00044 *> CMODE = 0 op2(C) = I 00045 *> CMODE = -1 op2(C) = inv(C) 00046 *> The Skeel condition number cond(A) = norminf( |inv(A)||A| ) 00047 *> is computed by computing scaling factors R such that 00048 *> diag(R)*A*op2(C) is row equilibrated and computing the standard 00049 *> infinity-norm condition number. 00050 *> \endverbatim 00051 * 00052 * Arguments: 00053 * ========== 00054 * 00055 *> \param[in] UPLO 00056 *> \verbatim 00057 *> UPLO is CHARACTER*1 00058 *> = 'U': Upper triangle of A is stored; 00059 *> = 'L': Lower triangle of A is stored. 00060 *> \endverbatim 00061 *> 00062 *> \param[in] N 00063 *> \verbatim 00064 *> N is INTEGER 00065 *> The number of linear equations, i.e., the order of the 00066 *> matrix A. N >= 0. 00067 *> \endverbatim 00068 *> 00069 *> \param[in] A 00070 *> \verbatim 00071 *> A is DOUBLE PRECISION array, dimension (LDA,N) 00072 *> On entry, the N-by-N matrix A. 00073 *> \endverbatim 00074 *> 00075 *> \param[in] LDA 00076 *> \verbatim 00077 *> LDA is INTEGER 00078 *> The leading dimension of the array A. LDA >= max(1,N). 00079 *> \endverbatim 00080 *> 00081 *> \param[in] AF 00082 *> \verbatim 00083 *> AF is DOUBLE PRECISION array, dimension (LDAF,N) 00084 *> The triangular factor U or L from the Cholesky factorization 00085 *> A = U**T*U or A = L*L**T, as computed by DPOTRF. 00086 *> \endverbatim 00087 *> 00088 *> \param[in] LDAF 00089 *> \verbatim 00090 *> LDAF is INTEGER 00091 *> The leading dimension of the array AF. LDAF >= max(1,N). 00092 *> \endverbatim 00093 *> 00094 *> \param[in] CMODE 00095 *> \verbatim 00096 *> CMODE is INTEGER 00097 *> Determines op2(C) in the formula op(A) * op2(C) as follows: 00098 *> CMODE = 1 op2(C) = C 00099 *> CMODE = 0 op2(C) = I 00100 *> CMODE = -1 op2(C) = inv(C) 00101 *> \endverbatim 00102 *> 00103 *> \param[in] C 00104 *> \verbatim 00105 *> C is DOUBLE PRECISION array, dimension (N) 00106 *> The vector C in the formula op(A) * op2(C). 00107 *> \endverbatim 00108 *> 00109 *> \param[out] INFO 00110 *> \verbatim 00111 *> INFO is INTEGER 00112 *> = 0: Successful exit. 00113 *> i > 0: The ith argument is invalid. 00114 *> \endverbatim 00115 *> 00116 *> \param[in] WORK 00117 *> \verbatim 00118 *> WORK is DOUBLE PRECISION array, dimension (3*N). 00119 *> Workspace. 00120 *> \endverbatim 00121 *> 00122 *> \param[in] IWORK 00123 *> \verbatim 00124 *> IWORK is INTEGER array, dimension (N). 00125 *> Workspace. 00126 *> \endverbatim 00127 * 00128 * Authors: 00129 * ======== 00130 * 00131 *> \author Univ. of Tennessee 00132 *> \author Univ. of California Berkeley 00133 *> \author Univ. of Colorado Denver 00134 *> \author NAG Ltd. 00135 * 00136 *> \date November 2011 00137 * 00138 *> \ingroup doublePOcomputational 00139 * 00140 * ===================================================================== 00141 DOUBLE PRECISION FUNCTION DLA_PORCOND( UPLO, N, A, LDA, AF, LDAF, 00142 $ CMODE, C, INFO, WORK, 00143 $ IWORK ) 00144 * 00145 * -- LAPACK computational routine (version 3.4.0) -- 00146 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00147 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00148 * November 2011 00149 * 00150 * .. Scalar Arguments .. 00151 CHARACTER UPLO 00152 INTEGER N, LDA, LDAF, INFO, CMODE 00153 DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), WORK( * ), 00154 $ C( * ) 00155 * .. 00156 * .. Array Arguments .. 00157 INTEGER IWORK( * ) 00158 * .. 00159 * 00160 * ===================================================================== 00161 * 00162 * .. Local Scalars .. 00163 INTEGER KASE, I, J 00164 DOUBLE PRECISION AINVNM, TMP 00165 LOGICAL UP 00166 * .. 00167 * .. Array Arguments .. 00168 INTEGER ISAVE( 3 ) 00169 * .. 00170 * .. External Functions .. 00171 LOGICAL LSAME 00172 INTEGER IDAMAX 00173 EXTERNAL LSAME, IDAMAX 00174 * .. 00175 * .. External Subroutines .. 00176 EXTERNAL DLACN2, DPOTRS, XERBLA 00177 * .. 00178 * .. Intrinsic Functions .. 00179 INTRINSIC ABS, MAX 00180 * .. 00181 * .. Executable Statements .. 00182 * 00183 DLA_PORCOND = 0.0D+0 00184 * 00185 INFO = 0 00186 IF( N.LT.0 ) THEN 00187 INFO = -2 00188 END IF 00189 IF( INFO.NE.0 ) THEN 00190 CALL XERBLA( 'DLA_PORCOND', -INFO ) 00191 RETURN 00192 END IF 00193 00194 IF( N.EQ.0 ) THEN 00195 DLA_PORCOND = 1.0D+0 00196 RETURN 00197 END IF 00198 UP = .FALSE. 00199 IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE. 00200 * 00201 * Compute the equilibration matrix R such that 00202 * inv(R)*A*C has unit 1-norm. 00203 * 00204 IF ( UP ) THEN 00205 DO I = 1, N 00206 TMP = 0.0D+0 00207 IF ( CMODE .EQ. 1 ) THEN 00208 DO J = 1, I 00209 TMP = TMP + ABS( A( J, I ) * C( J ) ) 00210 END DO 00211 DO J = I+1, N 00212 TMP = TMP + ABS( A( I, J ) * C( J ) ) 00213 END DO 00214 ELSE IF ( CMODE .EQ. 0 ) THEN 00215 DO J = 1, I 00216 TMP = TMP + ABS( A( J, I ) ) 00217 END DO 00218 DO J = I+1, N 00219 TMP = TMP + ABS( A( I, J ) ) 00220 END DO 00221 ELSE 00222 DO J = 1, I 00223 TMP = TMP + ABS( A( J ,I ) / C( J ) ) 00224 END DO 00225 DO J = I+1, N 00226 TMP = TMP + ABS( A( I, J ) / C( J ) ) 00227 END DO 00228 END IF 00229 WORK( 2*N+I ) = TMP 00230 END DO 00231 ELSE 00232 DO I = 1, N 00233 TMP = 0.0D+0 00234 IF ( CMODE .EQ. 1 ) THEN 00235 DO J = 1, I 00236 TMP = TMP + ABS( A( I, J ) * C( J ) ) 00237 END DO 00238 DO J = I+1, N 00239 TMP = TMP + ABS( A( J, I ) * C( J ) ) 00240 END DO 00241 ELSE IF ( CMODE .EQ. 0 ) THEN 00242 DO J = 1, I 00243 TMP = TMP + ABS( A( I, J ) ) 00244 END DO 00245 DO J = I+1, N 00246 TMP = TMP + ABS( A( J, I ) ) 00247 END DO 00248 ELSE 00249 DO J = 1, I 00250 TMP = TMP + ABS( A( I, J ) / C( J ) ) 00251 END DO 00252 DO J = I+1, N 00253 TMP = TMP + ABS( A( J, I ) / C( J ) ) 00254 END DO 00255 END IF 00256 WORK( 2*N+I ) = TMP 00257 END DO 00258 ENDIF 00259 * 00260 * Estimate the norm of inv(op(A)). 00261 * 00262 AINVNM = 0.0D+0 00263 00264 KASE = 0 00265 10 CONTINUE 00266 CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE ) 00267 IF( KASE.NE.0 ) THEN 00268 IF( KASE.EQ.2 ) THEN 00269 * 00270 * Multiply by R. 00271 * 00272 DO I = 1, N 00273 WORK( I ) = WORK( I ) * WORK( 2*N+I ) 00274 END DO 00275 00276 IF (UP) THEN 00277 CALL DPOTRS( 'Upper', N, 1, AF, LDAF, WORK, N, INFO ) 00278 ELSE 00279 CALL DPOTRS( 'Lower', N, 1, AF, LDAF, WORK, N, INFO ) 00280 ENDIF 00281 * 00282 * Multiply by inv(C). 00283 * 00284 IF ( CMODE .EQ. 1 ) THEN 00285 DO I = 1, N 00286 WORK( I ) = WORK( I ) / C( I ) 00287 END DO 00288 ELSE IF ( CMODE .EQ. -1 ) THEN 00289 DO I = 1, N 00290 WORK( I ) = WORK( I ) * C( I ) 00291 END DO 00292 END IF 00293 ELSE 00294 * 00295 * Multiply by inv(C**T). 00296 * 00297 IF ( CMODE .EQ. 1 ) THEN 00298 DO I = 1, N 00299 WORK( I ) = WORK( I ) / C( I ) 00300 END DO 00301 ELSE IF ( CMODE .EQ. -1 ) THEN 00302 DO I = 1, N 00303 WORK( I ) = WORK( I ) * C( I ) 00304 END DO 00305 END IF 00306 00307 IF ( UP ) THEN 00308 CALL DPOTRS( 'Upper', N, 1, AF, LDAF, WORK, N, INFO ) 00309 ELSE 00310 CALL DPOTRS( 'Lower', N, 1, AF, LDAF, WORK, N, INFO ) 00311 ENDIF 00312 * 00313 * Multiply by R. 00314 * 00315 DO I = 1, N 00316 WORK( I ) = WORK( I ) * WORK( 2*N+I ) 00317 END DO 00318 END IF 00319 GO TO 10 00320 END IF 00321 * 00322 * Compute the estimate of the reciprocal condition number. 00323 * 00324 IF( AINVNM .NE. 0.0D+0 ) 00325 $ DLA_PORCOND = ( 1.0D+0 / AINVNM ) 00326 * 00327 RETURN 00328 * 00329 END