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