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