![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zla_porcond_x.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_x.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: DOUBLE PRECISION FUNCTION ZLA_PORCOND_X( UPLO, N, A, LDA, AF,
2: $ LDAF, X, INFO, WORK,
3: $ 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: INTEGER N, LDA, LDAF, INFO
18: * ..
19: * .. Array Arguments ..
20: COMPLEX*16 A( LDA, * ), AF( LDAF, * ), WORK( * ), X( * )
21: DOUBLE PRECISION RWORK( * )
22: * ..
23: *
24: * Purpose
25: * =======
26: *
27: * ZLA_PORCOND_X Computes the infinity norm condition number of
28: * op(A) * diag(X) where X is a COMPLEX*16 vector.
29: *
30: * Arguments
31: * =========
32: *
33: * UPLO (input) CHARACTER*1
34: * = 'U': Upper triangle of A is stored;
35: * = 'L': Lower triangle of A is stored.
36: *
37: * N (input) INTEGER
38: * The number of linear equations, i.e., the order of the
39: * matrix A. N >= 0.
40: *
41: * A (input) COMPLEX*16 array, dimension (LDA,N)
42: * On entry, the N-by-N matrix A.
43: *
44: * LDA (input) INTEGER
45: * The leading dimension of the array A. LDA >= max(1,N).
46: *
47: * AF (input) COMPLEX*16 array, dimension (LDAF,N)
48: * The triangular factor U or L from the Cholesky factorization
49: * A = U**T*U or A = L*L**T, as computed by ZPOTRF.
50: *
51: * LDAF (input) INTEGER
52: * The leading dimension of the array AF. LDAF >= max(1,N).
53: *
54: * X (input) COMPLEX*16 array, dimension (N)
55: * The vector X in the formula op(A) * diag(X).
56: *
57: * INFO (output) INTEGER
58: * = 0: Successful exit.
59: * i > 0: The ith argument is invalid.
60: *
61: * WORK (input) COMPLEX*16 array, dimension (2*N).
62: * Workspace.
63: *
64: * RWORK (input) DOUBLE PRECISION array, dimension (N).
65: * Workspace.
66: *
67: * =====================================================================
68: *
69: * .. Local Scalars ..
70: INTEGER KASE, I, J
71: DOUBLE PRECISION AINVNM, ANORM, TMP
72: LOGICAL UP
73: COMPLEX*16 ZDUM
74: * ..
75: * .. Local Arrays ..
76: INTEGER ISAVE( 3 )
77: * ..
78: * .. External Functions ..
79: LOGICAL LSAME
80: EXTERNAL LSAME
81: * ..
82: * .. External Subroutines ..
83: EXTERNAL ZLACN2, ZPOTRS, XERBLA
84: * ..
85: * .. Intrinsic Functions ..
86: INTRINSIC ABS, MAX, REAL, DIMAG
87: * ..
88: * .. Statement Functions ..
89: DOUBLE PRECISION CABS1
90: * ..
91: * .. Statement Function Definitions ..
92: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
93: * ..
94: * .. Executable Statements ..
95: *
96: ZLA_PORCOND_X = 0.0D+0
97: *
98: INFO = 0
99: IF( N.LT.0 ) THEN
100: INFO = -2
101: END IF
102: IF( INFO.NE.0 ) THEN
103: CALL XERBLA( 'ZLA_PORCOND_X', -INFO )
104: RETURN
105: END IF
106: UP = .FALSE.
107: IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE.
108: *
109: * Compute norm of op(A)*op2(C).
110: *
111: ANORM = 0.0D+0
112: IF ( UP ) THEN
113: DO I = 1, N
114: TMP = 0.0D+0
115: DO J = 1, I
116: TMP = TMP + CABS1( A( J, I ) * X( J ) )
117: END DO
118: DO J = I+1, N
119: TMP = TMP + CABS1( A( I, J ) * X( J ) )
120: END DO
121: RWORK( I ) = TMP
122: ANORM = MAX( ANORM, TMP )
123: END DO
124: ELSE
125: DO I = 1, N
126: TMP = 0.0D+0
127: DO J = 1, I
128: TMP = TMP + CABS1( A( I, J ) * X( J ) )
129: END DO
130: DO J = I+1, N
131: TMP = TMP + CABS1( A( J, I ) * X( J ) )
132: END DO
133: RWORK( I ) = TMP
134: ANORM = MAX( ANORM, TMP )
135: END DO
136: END IF
137: *
138: * Quick return if possible.
139: *
140: IF( N.EQ.0 ) THEN
141: ZLA_PORCOND_X = 1.0D+0
142: RETURN
143: ELSE IF( ANORM .EQ. 0.0D+0 ) THEN
144: RETURN
145: END IF
146: *
147: * Estimate the norm of inv(op(A)).
148: *
149: AINVNM = 0.0D+0
150: *
151: KASE = 0
152: 10 CONTINUE
153: CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
154: IF( KASE.NE.0 ) THEN
155: IF( KASE.EQ.2 ) THEN
156: *
157: * Multiply by R.
158: *
159: DO I = 1, N
160: WORK( I ) = WORK( I ) * RWORK( I )
161: END DO
162: *
163: IF ( UP ) THEN
164: CALL ZPOTRS( 'U', N, 1, AF, LDAF,
165: $ WORK, N, INFO )
166: ELSE
167: CALL ZPOTRS( 'L', N, 1, AF, LDAF,
168: $ WORK, N, INFO )
169: ENDIF
170: *
171: * Multiply by inv(X).
172: *
173: DO I = 1, N
174: WORK( I ) = WORK( I ) / X( I )
175: END DO
176: ELSE
177: *
178: * Multiply by inv(X').
179: *
180: DO I = 1, N
181: WORK( I ) = WORK( I ) / X( I )
182: END DO
183: *
184: IF ( UP ) THEN
185: CALL ZPOTRS( 'U', N, 1, AF, LDAF,
186: $ WORK, N, INFO )
187: ELSE
188: CALL ZPOTRS( 'L', N, 1, AF, LDAF,
189: $ WORK, N, INFO )
190: END IF
191: *
192: * Multiply by R.
193: *
194: DO I = 1, N
195: WORK( I ) = WORK( I ) * RWORK( I )
196: END DO
197: END IF
198: GO TO 10
199: END IF
200: *
201: * Compute the estimate of the reciprocal condition number.
202: *
203: IF( AINVNM .NE. 0.0D+0 )
204: $ ZLA_PORCOND_X = 1.0D+0 / AINVNM
205: *
206: RETURN
207: *
208: END