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