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