![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zla_hercond_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_hercond_x.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: DOUBLE PRECISION FUNCTION ZLA_HERCOND_X( UPLO, N, A, LDA, AF,
2: $ LDAF, IPIV, X, 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: INTEGER N, LDA, LDAF, INFO
18: * ..
19: * .. Array Arguments ..
20: INTEGER IPIV( * )
21: COMPLEX*16 A( LDA, * ), AF( LDAF, * ), WORK( * ), X( * )
22: DOUBLE PRECISION RWORK( * )
23: * ..
24: *
25: * Purpose
26: * =======
27: *
28: * ZLA_HERCOND_X computes the infinity norm condition number of
29: * op(A) * diag(X) where X is a COMPLEX*16 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 block diagonal matrix D and the multipliers used to
50: * obtain the factor U or L as computed by ZHETRF.
51: *
52: * LDAF (input) INTEGER
53: * The leading dimension of the array AF. LDAF >= max(1,N).
54: *
55: * IPIV (input) INTEGER array, dimension (N)
56: * Details of the interchanges and the block structure of D
57: * as determined by CHETRF.
58: *
59: * X (input) COMPLEX*16 array, dimension (N)
60: * The vector X in the formula op(A) * diag(X).
61: *
62: * INFO (output) INTEGER
63: * = 0: Successful exit.
64: * i > 0: The ith argument is invalid.
65: *
66: * WORK (input) COMPLEX*16 array, dimension (2*N).
67: * Workspace.
68: *
69: * RWORK (input) DOUBLE PRECISION array, dimension (N).
70: * Workspace.
71: *
72: * =====================================================================
73: *
74: * .. Local Scalars ..
75: INTEGER KASE, I, J
76: DOUBLE PRECISION AINVNM, ANORM, TMP
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, ZHETRS, XERBLA
89: * ..
90: * .. Intrinsic Functions ..
91: INTRINSIC ABS, MAX
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_HERCOND_X = 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_HERCOND_X', -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: DO J = 1, I
121: TMP = TMP + CABS1( A( J, I ) * X( J ) )
122: END DO
123: DO J = I+1, N
124: TMP = TMP + CABS1( A( I, J ) * X( J ) )
125: END DO
126: RWORK( I ) = TMP
127: ANORM = MAX( ANORM, TMP )
128: END DO
129: ELSE
130: DO I = 1, N
131: TMP = 0.0D+0
132: DO J = 1, I
133: TMP = TMP + CABS1( A( I, J ) * X( J ) )
134: END DO
135: DO J = I+1, N
136: TMP = TMP + CABS1( A( J, I ) * X( J ) )
137: END DO
138: RWORK( I ) = TMP
139: ANORM = MAX( ANORM, TMP )
140: END DO
141: END IF
142: *
143: * Quick return if possible.
144: *
145: IF( N.EQ.0 ) THEN
146: ZLA_HERCOND_X = 1.0D+0
147: RETURN
148: ELSE IF( ANORM .EQ. 0.0D+0 ) THEN
149: RETURN
150: END IF
151: *
152: * Estimate the norm of inv(op(A)).
153: *
154: AINVNM = 0.0D+0
155: *
156: KASE = 0
157: 10 CONTINUE
158: CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
159: IF( KASE.NE.0 ) THEN
160: IF( KASE.EQ.2 ) THEN
161: *
162: * Multiply by R.
163: *
164: DO I = 1, N
165: WORK( I ) = WORK( I ) * RWORK( I )
166: END DO
167: *
168: IF ( UP ) THEN
169: CALL ZHETRS( 'U', N, 1, AF, LDAF, IPIV,
170: $ WORK, N, INFO )
171: ELSE
172: CALL ZHETRS( 'L', N, 1, AF, LDAF, IPIV,
173: $ WORK, N, INFO )
174: ENDIF
175: *
176: * Multiply by inv(X).
177: *
178: DO I = 1, N
179: WORK( I ) = WORK( I ) / X( I )
180: END DO
181: ELSE
182: *
183: * Multiply by inv(X').
184: *
185: DO I = 1, N
186: WORK( I ) = WORK( I ) / X( I )
187: END DO
188: *
189: IF ( UP ) THEN
190: CALL ZHETRS( 'U', N, 1, AF, LDAF, IPIV,
191: $ WORK, N, INFO )
192: ELSE
193: CALL ZHETRS( 'L', N, 1, AF, LDAF, IPIV,
194: $ WORK, N, INFO )
195: END IF
196: *
197: * Multiply by R.
198: *
199: DO I = 1, N
200: WORK( I ) = WORK( I ) * RWORK( I )
201: END DO
202: END IF
203: GO TO 10
204: END IF
205: *
206: * Compute the estimate of the reciprocal condition number.
207: *
208: IF( AINVNM .NE. 0.0D+0 )
209: $ ZLA_HERCOND_X = 1.0D+0 / AINVNM
210: *
211: RETURN
212: *
213: END