![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zla_porcond_c.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to zla_porcond_c.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: DOUBLE PRECISION FUNCTION ZLA_PORCOND_C( UPLO, N, A, LDA, AF,
2: $ LDAF, C, CAPPLY, INFO,
3: $ WORK, RWORK )
4: *
5: * -- LAPACK routine (version 3.2.1) --
6: * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
7: * -- Jason Riedy of Univ. of California Berkeley. --
8: * -- April 2009 --
9: *
10: * -- LAPACK is a software package provided by Univ. of Tennessee, --
11: * -- Univ. of California Berkeley and NAG Ltd. --
12: *
13: IMPLICIT NONE
14: * ..
15: * .. Scalar Arguments ..
16: CHARACTER UPLO
17: LOGICAL CAPPLY
18: INTEGER N, LDA, LDAF, INFO
19: * ..
20: * .. Array Arguments ..
21: COMPLEX*16 A( LDA, * ), AF( LDAF, * ), WORK( * )
22: DOUBLE PRECISION C( * ), RWORK( * )
23: * ..
24: *
25: * Purpose
26: * =======
27: *
28: * ZLA_PORCOND_C Computes the infinity norm condition number of
29: * op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector
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 number of linear equations, i.e., the order of the
40: * matrix A. N >= 0.
41: *
42: * A (input) COMPLEX*16 array, dimension (LDA,N)
43: * On entry, the N-by-N matrix A
44: *
45: * LDA (input) INTEGER
46: * The leading dimension of the array A. LDA >= max(1,N).
47: *
48: * AF (input) COMPLEX*16 array, dimension (LDAF,N)
49: * The triangular factor U or L from the Cholesky factorization
1.5 ! bertrand 50: * A = U**H*U or A = L*L**H, as computed by ZPOTRF.
1.1 bertrand 51: *
52: * LDAF (input) INTEGER
53: * The leading dimension of the array AF. LDAF >= max(1,N).
54: *
55: * C (input) DOUBLE PRECISION array, dimension (N)
56: * The vector C in the formula op(A) * inv(diag(C)).
57: *
58: * CAPPLY (input) LOGICAL
59: * If .TRUE. then access the vector C in the formula above.
60: *
61: * INFO (output) INTEGER
62: * = 0: Successful exit.
63: * i > 0: The ith argument is invalid.
64: *
65: * WORK (input) COMPLEX*16 array, dimension (2*N).
66: * Workspace.
67: *
68: * RWORK (input) DOUBLE PRECISION array, dimension (N).
69: * Workspace.
70: *
71: * =====================================================================
72: *
73: * .. Local Scalars ..
74: INTEGER KASE
75: DOUBLE PRECISION AINVNM, ANORM, TMP
76: INTEGER I, J
77: LOGICAL UP
78: COMPLEX*16 ZDUM
79: * ..
80: * .. Local Arrays ..
81: INTEGER ISAVE( 3 )
82: * ..
83: * .. External Functions ..
84: LOGICAL LSAME
85: EXTERNAL LSAME
86: * ..
87: * .. External Subroutines ..
88: EXTERNAL ZLACN2, ZPOTRS, XERBLA
89: * ..
90: * .. Intrinsic Functions ..
91: INTRINSIC ABS, MAX, REAL, 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: ZLA_PORCOND_C = 0.0D+0
102: *
103: INFO = 0
104: IF( N.LT.0 ) THEN
105: INFO = -2
106: END IF
107: IF( INFO.NE.0 ) THEN
108: CALL XERBLA( 'ZLA_PORCOND_C', -INFO )
109: RETURN
110: END IF
111: UP = .FALSE.
112: IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE.
113: *
114: * Compute norm of op(A)*op2(C).
115: *
116: ANORM = 0.0D+0
117: IF ( UP ) THEN
118: DO I = 1, N
119: TMP = 0.0D+0
120: IF ( CAPPLY ) THEN
121: DO J = 1, I
122: TMP = TMP + CABS1( A( J, I ) ) / C( J )
123: END DO
124: DO J = I+1, N
125: TMP = TMP + CABS1( A( I, J ) ) / C( J )
126: END DO
127: ELSE
128: DO J = 1, I
129: TMP = TMP + CABS1( A( J, I ) )
130: END DO
131: DO J = I+1, N
132: TMP = TMP + CABS1( A( I, J ) )
133: END DO
134: END IF
135: RWORK( I ) = TMP
136: ANORM = MAX( ANORM, TMP )
137: END DO
138: ELSE
139: DO I = 1, N
140: TMP = 0.0D+0
141: IF ( CAPPLY ) THEN
142: DO J = 1, I
143: TMP = TMP + CABS1( A( I, J ) ) / C( J )
144: END DO
145: DO J = I+1, N
146: TMP = TMP + CABS1( A( J, I ) ) / C( J )
147: END DO
148: ELSE
149: DO J = 1, I
150: TMP = TMP + CABS1( A( I, J ) )
151: END DO
152: DO J = I+1, N
153: TMP = TMP + CABS1( A( J, I ) )
154: END DO
155: END IF
156: RWORK( I ) = TMP
157: ANORM = MAX( ANORM, TMP )
158: END DO
159: END IF
160: *
161: * Quick return if possible.
162: *
163: IF( N.EQ.0 ) THEN
164: ZLA_PORCOND_C = 1.0D+0
165: RETURN
166: ELSE IF( ANORM .EQ. 0.0D+0 ) THEN
167: RETURN
168: END IF
169: *
170: * Estimate the norm of inv(op(A)).
171: *
172: AINVNM = 0.0D+0
173: *
174: KASE = 0
175: 10 CONTINUE
176: CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
177: IF( KASE.NE.0 ) THEN
178: IF( KASE.EQ.2 ) THEN
179: *
180: * Multiply by R.
181: *
182: DO I = 1, N
183: WORK( I ) = WORK( I ) * RWORK( I )
184: END DO
185: *
186: IF ( UP ) THEN
187: CALL ZPOTRS( 'U', N, 1, AF, LDAF,
188: $ WORK, N, INFO )
189: ELSE
190: CALL ZPOTRS( 'L', N, 1, AF, LDAF,
191: $ WORK, N, INFO )
192: ENDIF
193: *
194: * Multiply by inv(C).
195: *
196: IF ( CAPPLY ) THEN
197: DO I = 1, N
198: WORK( I ) = WORK( I ) * C( I )
199: END DO
200: END IF
201: ELSE
202: *
1.5 ! bertrand 203: * Multiply by inv(C**H).
1.1 bertrand 204: *
205: IF ( CAPPLY ) THEN
206: DO I = 1, N
207: WORK( I ) = WORK( I ) * C( I )
208: END DO
209: END IF
210: *
211: IF ( UP ) THEN
212: CALL ZPOTRS( 'U', N, 1, AF, LDAF,
213: $ WORK, N, INFO )
214: ELSE
215: CALL ZPOTRS( 'L', N, 1, AF, LDAF,
216: $ WORK, N, INFO )
217: END IF
218: *
219: * Multiply by R.
220: *
221: DO I = 1, N
222: WORK( I ) = WORK( I ) * RWORK( I )
223: END DO
224: END IF
225: GO TO 10
226: END IF
227: *
228: * Compute the estimate of the reciprocal condition number.
229: *
230: IF( AINVNM .NE. 0.0D+0 )
231: $ ZLA_PORCOND_C = 1.0D+0 / AINVNM
232: *
233: RETURN
234: *
235: END