1*> \brief \b DLA_GERCOND estimates the Skeel condition number for a general matrix. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DLA_GERCOND + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dla_gercond.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dla_gercond.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dla_gercond.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* DOUBLE PRECISION FUNCTION DLA_GERCOND( TRANS, N, A, LDA, AF, 22* LDAF, IPIV, CMODE, C, 23* INFO, WORK, IWORK ) 24* 25* .. Scalar Arguments .. 26* CHARACTER TRANS 27* INTEGER N, LDA, LDAF, INFO, CMODE 28* .. 29* .. Array Arguments .. 30* INTEGER IPIV( * ), IWORK( * ) 31* DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), WORK( * ), 32* $ C( * ) 33* .. 34* 35* 36*> \par Purpose: 37* ============= 38*> 39*> \verbatim 40*> 41*> DLA_GERCOND estimates the Skeel condition number of op(A) * op2(C) 42*> where op2 is determined by CMODE as follows 43*> CMODE = 1 op2(C) = C 44*> CMODE = 0 op2(C) = I 45*> CMODE = -1 op2(C) = inv(C) 46*> The Skeel condition number cond(A) = norminf( |inv(A)||A| ) 47*> is computed by computing scaling factors R such that 48*> diag(R)*A*op2(C) is row equilibrated and computing the standard 49*> infinity-norm condition number. 50*> \endverbatim 51* 52* Arguments: 53* ========== 54* 55*> \param[in] TRANS 56*> \verbatim 57*> TRANS is CHARACTER*1 58*> Specifies the form of the system of equations: 59*> = 'N': A * X = B (No transpose) 60*> = 'T': A**T * X = B (Transpose) 61*> = 'C': A**H * X = B (Conjugate Transpose = Transpose) 62*> \endverbatim 63*> 64*> \param[in] N 65*> \verbatim 66*> N is INTEGER 67*> The number of linear equations, i.e., the order of the 68*> matrix A. N >= 0. 69*> \endverbatim 70*> 71*> \param[in] A 72*> \verbatim 73*> A is DOUBLE PRECISION array, dimension (LDA,N) 74*> On entry, the N-by-N matrix A. 75*> \endverbatim 76*> 77*> \param[in] LDA 78*> \verbatim 79*> LDA is INTEGER 80*> The leading dimension of the array A. LDA >= max(1,N). 81*> \endverbatim 82*> 83*> \param[in] AF 84*> \verbatim 85*> AF is DOUBLE PRECISION array, dimension (LDAF,N) 86*> The factors L and U from the factorization 87*> A = P*L*U as computed by DGETRF. 88*> \endverbatim 89*> 90*> \param[in] LDAF 91*> \verbatim 92*> LDAF is INTEGER 93*> The leading dimension of the array AF. LDAF >= max(1,N). 94*> \endverbatim 95*> 96*> \param[in] IPIV 97*> \verbatim 98*> IPIV is INTEGER array, dimension (N) 99*> The pivot indices from the factorization A = P*L*U 100*> as computed by DGETRF; row i of the matrix was interchanged 101*> with row IPIV(i). 102*> \endverbatim 103*> 104*> \param[in] CMODE 105*> \verbatim 106*> CMODE is INTEGER 107*> Determines op2(C) in the formula op(A) * op2(C) as follows: 108*> CMODE = 1 op2(C) = C 109*> CMODE = 0 op2(C) = I 110*> CMODE = -1 op2(C) = inv(C) 111*> \endverbatim 112*> 113*> \param[in] C 114*> \verbatim 115*> C is DOUBLE PRECISION array, dimension (N) 116*> The vector C in the formula op(A) * op2(C). 117*> \endverbatim 118*> 119*> \param[out] INFO 120*> \verbatim 121*> INFO is INTEGER 122*> = 0: Successful exit. 123*> i > 0: The ith argument is invalid. 124*> \endverbatim 125*> 126*> \param[out] WORK 127*> \verbatim 128*> WORK is DOUBLE PRECISION array, dimension (3*N). 129*> Workspace. 130*> \endverbatim 131*> 132*> \param[out] IWORK 133*> \verbatim 134*> IWORK is INTEGER array, dimension (N). 135*> Workspace. 136*> \endverbatim 137* 138* Authors: 139* ======== 140* 141*> \author Univ. of Tennessee 142*> \author Univ. of California Berkeley 143*> \author Univ. of Colorado Denver 144*> \author NAG Ltd. 145* 146*> \ingroup doubleGEcomputational 147* 148* ===================================================================== 149 DOUBLE PRECISION FUNCTION DLA_GERCOND( TRANS, N, A, LDA, AF, 150 $ LDAF, IPIV, CMODE, C, 151 $ INFO, WORK, IWORK ) 152* 153* -- LAPACK computational routine -- 154* -- LAPACK is a software package provided by Univ. of Tennessee, -- 155* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 156* 157* .. Scalar Arguments .. 158 CHARACTER TRANS 159 INTEGER N, LDA, LDAF, INFO, CMODE 160* .. 161* .. Array Arguments .. 162 INTEGER IPIV( * ), IWORK( * ) 163 DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), WORK( * ), 164 $ C( * ) 165* .. 166* 167* ===================================================================== 168* 169* .. Local Scalars .. 170 LOGICAL NOTRANS 171 INTEGER KASE, I, J 172 DOUBLE PRECISION AINVNM, TMP 173* .. 174* .. Local Arrays .. 175 INTEGER ISAVE( 3 ) 176* .. 177* .. External Functions .. 178 LOGICAL LSAME 179 EXTERNAL LSAME 180* .. 181* .. External Subroutines .. 182 EXTERNAL DLACN2, DGETRS, XERBLA 183* .. 184* .. Intrinsic Functions .. 185 INTRINSIC ABS, MAX 186* .. 187* .. Executable Statements .. 188* 189 DLA_GERCOND = 0.0D+0 190* 191 INFO = 0 192 NOTRANS = LSAME( TRANS, 'N' ) 193 IF ( .NOT. NOTRANS .AND. .NOT. LSAME(TRANS, 'T') 194 $ .AND. .NOT. LSAME(TRANS, 'C') ) THEN 195 INFO = -1 196 ELSE IF( N.LT.0 ) THEN 197 INFO = -2 198 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 199 INFO = -4 200 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN 201 INFO = -6 202 END IF 203 IF( INFO.NE.0 ) THEN 204 CALL XERBLA( 'DLA_GERCOND', -INFO ) 205 RETURN 206 END IF 207 IF( N.EQ.0 ) THEN 208 DLA_GERCOND = 1.0D+0 209 RETURN 210 END IF 211* 212* Compute the equilibration matrix R such that 213* inv(R)*A*C has unit 1-norm. 214* 215 IF (NOTRANS) THEN 216 DO I = 1, N 217 TMP = 0.0D+0 218 IF ( CMODE .EQ. 1 ) THEN 219 DO J = 1, N 220 TMP = TMP + ABS( A( I, J ) * C( J ) ) 221 END DO 222 ELSE IF ( CMODE .EQ. 0 ) THEN 223 DO J = 1, N 224 TMP = TMP + ABS( A( I, J ) ) 225 END DO 226 ELSE 227 DO J = 1, N 228 TMP = TMP + ABS( A( I, J ) / C( J ) ) 229 END DO 230 END IF 231 WORK( 2*N+I ) = TMP 232 END DO 233 ELSE 234 DO I = 1, N 235 TMP = 0.0D+0 236 IF ( CMODE .EQ. 1 ) THEN 237 DO J = 1, N 238 TMP = TMP + ABS( A( J, I ) * C( J ) ) 239 END DO 240 ELSE IF ( CMODE .EQ. 0 ) THEN 241 DO J = 1, N 242 TMP = TMP + ABS( A( J, I ) ) 243 END DO 244 ELSE 245 DO J = 1, N 246 TMP = TMP + ABS( A( J, I ) / C( J ) ) 247 END DO 248 END IF 249 WORK( 2*N+I ) = TMP 250 END DO 251 END IF 252* 253* Estimate the norm of inv(op(A)). 254* 255 AINVNM = 0.0D+0 256 257 KASE = 0 258 10 CONTINUE 259 CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE ) 260 IF( KASE.NE.0 ) THEN 261 IF( KASE.EQ.2 ) THEN 262* 263* Multiply by R. 264* 265 DO I = 1, N 266 WORK(I) = WORK(I) * WORK(2*N+I) 267 END DO 268 269 IF (NOTRANS) THEN 270 CALL DGETRS( 'No transpose', N, 1, AF, LDAF, IPIV, 271 $ WORK, N, INFO ) 272 ELSE 273 CALL DGETRS( 'Transpose', N, 1, AF, LDAF, IPIV, 274 $ WORK, N, INFO ) 275 END IF 276* 277* Multiply by inv(C). 278* 279 IF ( CMODE .EQ. 1 ) THEN 280 DO I = 1, N 281 WORK( I ) = WORK( I ) / C( I ) 282 END DO 283 ELSE IF ( CMODE .EQ. -1 ) THEN 284 DO I = 1, N 285 WORK( I ) = WORK( I ) * C( I ) 286 END DO 287 END IF 288 ELSE 289* 290* Multiply by inv(C**T). 291* 292 IF ( CMODE .EQ. 1 ) THEN 293 DO I = 1, N 294 WORK( I ) = WORK( I ) / C( I ) 295 END DO 296 ELSE IF ( CMODE .EQ. -1 ) THEN 297 DO I = 1, N 298 WORK( I ) = WORK( I ) * C( I ) 299 END DO 300 END IF 301 302 IF (NOTRANS) THEN 303 CALL DGETRS( 'Transpose', N, 1, AF, LDAF, IPIV, 304 $ WORK, N, INFO ) 305 ELSE 306 CALL DGETRS( 'No transpose', N, 1, AF, LDAF, IPIV, 307 $ WORK, N, INFO ) 308 END IF 309* 310* Multiply by R. 311* 312 DO I = 1, N 313 WORK( I ) = WORK( I ) * WORK( 2*N+I ) 314 END DO 315 END IF 316 GO TO 10 317 END IF 318* 319* Compute the estimate of the reciprocal condition number. 320* 321 IF( AINVNM .NE. 0.0D+0 ) 322 $ DLA_GERCOND = ( 1.0D+0 / AINVNM ) 323* 324 RETURN 325* 326* End of DLA_GERCOND 327* 328 END 329