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