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