Annotation of rpl/lapack/lapack/dpocon.f, revision 1.4
1.1 bertrand 1: SUBROUTINE DPOCON( UPLO, N, A, LDA, ANORM, RCOND, WORK, IWORK,
2: $ INFO )
3: *
4: * -- LAPACK routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH.
10: *
11: * .. Scalar Arguments ..
12: CHARACTER UPLO
13: INTEGER INFO, LDA, N
14: DOUBLE PRECISION ANORM, RCOND
15: * ..
16: * .. Array Arguments ..
17: INTEGER IWORK( * )
18: DOUBLE PRECISION A( LDA, * ), WORK( * )
19: * ..
20: *
21: * Purpose
22: * =======
23: *
24: * DPOCON estimates the reciprocal of the condition number (in the
25: * 1-norm) of a real symmetric positive definite matrix using the
26: * Cholesky factorization A = U**T*U or A = L*L**T computed by DPOTRF.
27: *
28: * An estimate is obtained for norm(inv(A)), and the reciprocal of the
29: * condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
30: *
31: * Arguments
32: * =========
33: *
34: * UPLO (input) CHARACTER*1
35: * = 'U': Upper triangle of A is stored;
36: * = 'L': Lower triangle of A is stored.
37: *
38: * N (input) INTEGER
39: * The order of the matrix A. N >= 0.
40: *
41: * A (input) DOUBLE PRECISION array, dimension (LDA,N)
42: * The triangular factor U or L from the Cholesky factorization
43: * A = U**T*U or A = L*L**T, as computed by DPOTRF.
44: *
45: * LDA (input) INTEGER
46: * The leading dimension of the array A. LDA >= max(1,N).
47: *
48: * ANORM (input) DOUBLE PRECISION
49: * The 1-norm (or infinity-norm) of the symmetric matrix A.
50: *
51: * RCOND (output) DOUBLE PRECISION
52: * The reciprocal of the condition number of the matrix A,
53: * computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
54: * estimate of the 1-norm of inv(A) computed in this routine.
55: *
56: * WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
57: *
58: * IWORK (workspace) INTEGER array, dimension (N)
59: *
60: * INFO (output) INTEGER
61: * = 0: successful exit
62: * < 0: if INFO = -i, the i-th argument had an illegal value
63: *
64: * =====================================================================
65: *
66: * .. Parameters ..
67: DOUBLE PRECISION ONE, ZERO
68: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
69: * ..
70: * .. Local Scalars ..
71: LOGICAL UPPER
72: CHARACTER NORMIN
73: INTEGER IX, KASE
74: DOUBLE PRECISION AINVNM, SCALE, SCALEL, SCALEU, SMLNUM
75: * ..
76: * .. Local Arrays ..
77: INTEGER ISAVE( 3 )
78: * ..
79: * .. External Functions ..
80: LOGICAL LSAME
81: INTEGER IDAMAX
82: DOUBLE PRECISION DLAMCH
83: EXTERNAL LSAME, IDAMAX, DLAMCH
84: * ..
85: * .. External Subroutines ..
86: EXTERNAL DLACN2, DLATRS, DRSCL, XERBLA
87: * ..
88: * .. Intrinsic Functions ..
89: INTRINSIC ABS, MAX
90: * ..
91: * .. Executable Statements ..
92: *
93: * Test the input parameters.
94: *
95: INFO = 0
96: UPPER = LSAME( UPLO, 'U' )
97: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
98: INFO = -1
99: ELSE IF( N.LT.0 ) THEN
100: INFO = -2
101: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
102: INFO = -4
103: ELSE IF( ANORM.LT.ZERO ) THEN
104: INFO = -5
105: END IF
106: IF( INFO.NE.0 ) THEN
107: CALL XERBLA( 'DPOCON', -INFO )
108: RETURN
109: END IF
110: *
111: * Quick return if possible
112: *
113: RCOND = ZERO
114: IF( N.EQ.0 ) THEN
115: RCOND = ONE
116: RETURN
117: ELSE IF( ANORM.EQ.ZERO ) THEN
118: RETURN
119: END IF
120: *
121: SMLNUM = DLAMCH( 'Safe minimum' )
122: *
123: * Estimate the 1-norm of inv(A).
124: *
125: KASE = 0
126: NORMIN = 'N'
127: 10 CONTINUE
128: CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
129: IF( KASE.NE.0 ) THEN
130: IF( UPPER ) THEN
131: *
132: * Multiply by inv(U').
133: *
134: CALL DLATRS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N, A,
135: $ LDA, WORK, SCALEL, WORK( 2*N+1 ), INFO )
136: NORMIN = 'Y'
137: *
138: * Multiply by inv(U).
139: *
140: CALL DLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
141: $ A, LDA, WORK, SCALEU, WORK( 2*N+1 ), INFO )
142: ELSE
143: *
144: * Multiply by inv(L).
145: *
146: CALL DLATRS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N,
147: $ A, LDA, WORK, SCALEL, WORK( 2*N+1 ), INFO )
148: NORMIN = 'Y'
149: *
150: * Multiply by inv(L').
151: *
152: CALL DLATRS( 'Lower', 'Transpose', 'Non-unit', NORMIN, N, A,
153: $ LDA, WORK, SCALEU, WORK( 2*N+1 ), INFO )
154: END IF
155: *
156: * Multiply by 1/SCALE if doing so will not cause overflow.
157: *
158: SCALE = SCALEL*SCALEU
159: IF( SCALE.NE.ONE ) THEN
160: IX = IDAMAX( N, WORK, 1 )
161: IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
162: $ GO TO 20
163: CALL DRSCL( N, SCALE, WORK, 1 )
164: END IF
165: GO TO 10
166: END IF
167: *
168: * Compute the estimate of the reciprocal condition number.
169: *
170: IF( AINVNM.NE.ZERO )
171: $ RCOND = ( ONE / AINVNM ) / ANORM
172: *
173: 20 CONTINUE
174: RETURN
175: *
176: * End of DPOCON
177: *
178: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dpocon.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack