Annotation of rpl/lapack/lapack/zpocon.f, revision 1.5
1.1 bertrand 1: SUBROUTINE ZPOCON( UPLO, N, A, LDA, ANORM, RCOND, WORK, RWORK,
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 ZLACN2 in place of ZLACON, 10 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: DOUBLE PRECISION RWORK( * )
18: COMPLEX*16 A( LDA, * ), WORK( * )
19: * ..
20: *
21: * Purpose
22: * =======
23: *
24: * ZPOCON estimates the reciprocal of the condition number (in the
25: * 1-norm) of a complex Hermitian positive definite matrix using the
26: * Cholesky factorization A = U**H*U or A = L*L**H computed by ZPOTRF.
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) COMPLEX*16 array, dimension (LDA,N)
42: * The triangular factor U or L from the Cholesky factorization
43: * A = U**H*U or A = L*L**H, as computed by ZPOTRF.
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 Hermitian 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) COMPLEX*16 array, dimension (2*N)
57: *
58: * RWORK (workspace) DOUBLE PRECISION 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: COMPLEX*16 ZDUM
76: * ..
77: * .. Local Arrays ..
78: INTEGER ISAVE( 3 )
79: * ..
80: * .. External Functions ..
81: LOGICAL LSAME
82: INTEGER IZAMAX
83: DOUBLE PRECISION DLAMCH
84: EXTERNAL LSAME, IZAMAX, DLAMCH
85: * ..
86: * .. External Subroutines ..
87: EXTERNAL XERBLA, ZDRSCL, ZLACN2, ZLATRS
88: * ..
89: * .. Intrinsic Functions ..
90: INTRINSIC ABS, DBLE, DIMAG, MAX
91: * ..
92: * .. Statement Functions ..
93: DOUBLE PRECISION CABS1
94: * ..
95: * .. Statement Function definitions ..
96: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
97: * ..
98: * .. Executable Statements ..
99: *
100: * Test the input parameters.
101: *
102: INFO = 0
103: UPPER = LSAME( UPLO, 'U' )
104: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
105: INFO = -1
106: ELSE IF( N.LT.0 ) THEN
107: INFO = -2
108: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
109: INFO = -4
110: ELSE IF( ANORM.LT.ZERO ) THEN
111: INFO = -5
112: END IF
113: IF( INFO.NE.0 ) THEN
114: CALL XERBLA( 'ZPOCON', -INFO )
115: RETURN
116: END IF
117: *
118: * Quick return if possible
119: *
120: RCOND = ZERO
121: IF( N.EQ.0 ) THEN
122: RCOND = ONE
123: RETURN
124: ELSE IF( ANORM.EQ.ZERO ) THEN
125: RETURN
126: END IF
127: *
128: SMLNUM = DLAMCH( 'Safe minimum' )
129: *
130: * Estimate the 1-norm of inv(A).
131: *
132: KASE = 0
133: NORMIN = 'N'
134: 10 CONTINUE
135: CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
136: IF( KASE.NE.0 ) THEN
137: IF( UPPER ) THEN
138: *
139: * Multiply by inv(U').
140: *
141: CALL ZLATRS( 'Upper', 'Conjugate transpose', 'Non-unit',
142: $ NORMIN, N, A, LDA, WORK, SCALEL, RWORK, INFO )
143: NORMIN = 'Y'
144: *
145: * Multiply by inv(U).
146: *
147: CALL ZLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
148: $ A, LDA, WORK, SCALEU, RWORK, INFO )
149: ELSE
150: *
151: * Multiply by inv(L).
152: *
153: CALL ZLATRS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N,
154: $ A, LDA, WORK, SCALEL, RWORK, INFO )
155: NORMIN = 'Y'
156: *
157: * Multiply by inv(L').
158: *
159: CALL ZLATRS( 'Lower', 'Conjugate transpose', 'Non-unit',
160: $ NORMIN, N, A, LDA, WORK, SCALEU, RWORK, INFO )
161: END IF
162: *
163: * Multiply by 1/SCALE if doing so will not cause overflow.
164: *
165: SCALE = SCALEL*SCALEU
166: IF( SCALE.NE.ONE ) THEN
167: IX = IZAMAX( N, WORK, 1 )
168: IF( SCALE.LT.CABS1( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
169: $ GO TO 20
170: CALL ZDRSCL( N, SCALE, WORK, 1 )
171: END IF
172: GO TO 10
173: END IF
174: *
175: * Compute the estimate of the reciprocal condition number.
176: *
177: IF( AINVNM.NE.ZERO )
178: $ RCOND = ( ONE / AINVNM ) / ANORM
179: *
180: 20 CONTINUE
181: RETURN
182: *
183: * End of ZPOCON
184: *
185: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zpocon.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack