1 SUBROUTINE SPOCON( UPLO, N, A, LDA, ANORM, RCOND, WORK, IWORK, 2 $ INFO ) 3* 4* -- LAPACK routine (version 3.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* March 31, 1993 8* 9* .. Scalar Arguments .. 10 CHARACTER UPLO 11 INTEGER INFO, LDA, N 12 REAL ANORM, RCOND 13* .. 14* .. Array Arguments .. 15 INTEGER IWORK( * ) 16 REAL A( LDA, * ), WORK( * ) 17* .. 18* 19* Purpose 20* ======= 21* 22* SPOCON estimates the reciprocal of the condition number (in the 23* 1-norm) of a real symmetric positive definite matrix using the 24* Cholesky factorization A = U**T*U or A = L*L**T computed by SPOTRF. 25* 26* An estimate is obtained for norm(inv(A)), and the reciprocal of the 27* condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))). 28* 29* Arguments 30* ========= 31* 32* UPLO (input) CHARACTER*1 33* = 'U': Upper triangle of A is stored; 34* = 'L': Lower triangle of A is stored. 35* 36* N (input) INTEGER 37* The order of the matrix A. N >= 0. 38* 39* A (input) REAL array, dimension (LDA,N) 40* The triangular factor U or L from the Cholesky factorization 41* A = U**T*U or A = L*L**T, as computed by SPOTRF. 42* 43* LDA (input) INTEGER 44* The leading dimension of the array A. LDA >= max(1,N). 45* 46* ANORM (input) REAL 47* The 1-norm (or infinity-norm) of the symmetric matrix A. 48* 49* RCOND (output) REAL 50* The reciprocal of the condition number of the matrix A, 51* computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an 52* estimate of the 1-norm of inv(A) computed in this routine. 53* 54* WORK (workspace) REAL array, dimension (3*N) 55* 56* IWORK (workspace) INTEGER array, dimension (N) 57* 58* INFO (output) INTEGER 59* = 0: successful exit 60* < 0: if INFO = -i, the i-th argument had an illegal value 61* 62* ===================================================================== 63* 64* .. Parameters .. 65 REAL ONE, ZERO 66 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 67* .. 68* .. Local Scalars .. 69 LOGICAL UPPER 70 CHARACTER NORMIN 71 INTEGER IX, KASE 72 REAL AINVNM, SCALE, SCALEL, SCALEU, SMLNUM 73* .. 74* .. External Functions .. 75 LOGICAL LSAME 76 INTEGER ISAMAX 77 REAL SLAMCH 78 EXTERNAL LSAME, ISAMAX, SLAMCH 79* .. 80* .. External Subroutines .. 81 EXTERNAL SLACON, SLATRS, SRSCL, XERBLA 82* .. 83* .. Intrinsic Functions .. 84 INTRINSIC ABS, MAX 85* .. 86* .. Executable Statements .. 87* 88* Test the input parameters. 89* 90 INFO = 0 91 UPPER = LSAME( UPLO, 'U' ) 92 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 93 INFO = -1 94 ELSE IF( N.LT.0 ) THEN 95 INFO = -2 96 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 97 INFO = -4 98 ELSE IF( ANORM.LT.ZERO ) THEN 99 INFO = -5 100 END IF 101 IF( INFO.NE.0 ) THEN 102 CALL XERBLA( 'SPOCON', -INFO ) 103 RETURN 104 END IF 105* 106* Quick return if possible 107* 108 RCOND = ZERO 109 IF( N.EQ.0 ) THEN 110 RCOND = ONE 111 RETURN 112 ELSE IF( ANORM.EQ.ZERO ) THEN 113 RETURN 114 END IF 115* 116 SMLNUM = SLAMCH( 'Safe minimum' ) 117* 118* Estimate the 1-norm of inv(A). 119* 120 KASE = 0 121 NORMIN = 'N' 122 10 CONTINUE 123 CALL SLACON( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE ) 124 IF( KASE.NE.0 ) THEN 125 IF( UPPER ) THEN 126* 127* Multiply by inv(U'). 128* 129 CALL SLATRS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N, A, 130 $ LDA, WORK, SCALEL, WORK( 2*N+1 ), INFO ) 131 NORMIN = 'Y' 132* 133* Multiply by inv(U). 134* 135 CALL SLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N, 136 $ A, LDA, WORK, SCALEU, WORK( 2*N+1 ), INFO ) 137 ELSE 138* 139* Multiply by inv(L). 140* 141 CALL SLATRS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N, 142 $ A, LDA, WORK, SCALEL, WORK( 2*N+1 ), INFO ) 143 NORMIN = 'Y' 144* 145* Multiply by inv(L'). 146* 147 CALL SLATRS( 'Lower', 'Transpose', 'Non-unit', NORMIN, N, A, 148 $ LDA, WORK, SCALEU, WORK( 2*N+1 ), INFO ) 149 END IF 150* 151* Multiply by 1/SCALE if doing so will not cause overflow. 152* 153 SCALE = SCALEL*SCALEU 154 IF( SCALE.NE.ONE ) THEN 155 IX = ISAMAX( N, WORK, 1 ) 156 IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO ) 157 $ GO TO 20 158 CALL SRSCL( N, SCALE, WORK, 1 ) 159 END IF 160 GO TO 10 161 END IF 162* 163* Compute the estimate of the reciprocal condition number. 164* 165 IF( AINVNM.NE.ZERO ) 166 $ RCOND = ( ONE / AINVNM ) / ANORM 167* 168 20 CONTINUE 169 RETURN 170* 171* End of SPOCON 172* 173 END 174