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