![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zla_syrcond_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_syrcond_x.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: DOUBLE PRECISION FUNCTION ZLA_SYRCOND_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_SYRCOND_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 ZSYTRF.
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 ZSYTRF.
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
76: DOUBLE PRECISION AINVNM, ANORM, TMP
77: INTEGER I, J
78: LOGICAL UP
79: COMPLEX*16 ZDUM
80: * ..
81: * .. Local Arrays ..
82: INTEGER ISAVE( 3 )
83: * ..
84: * .. External Functions ..
85: LOGICAL LSAME
86: EXTERNAL LSAME
87: * ..
88: * .. External Subroutines ..
89: EXTERNAL ZLACN2, ZSYTRS, XERBLA
90: * ..
91: * .. Intrinsic Functions ..
92: INTRINSIC ABS, MAX
93: * ..
94: * .. Statement Functions ..
95: DOUBLE PRECISION CABS1
96: * ..
97: * .. Statement Function Definitions ..
98: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
99: * ..
100: * .. Executable Statements ..
101: *
102: ZLA_SYRCOND_X = 0.0D+0
103: *
104: INFO = 0
105: IF( N.LT.0 ) THEN
106: INFO = -2
107: END IF
108: IF( INFO.NE.0 ) THEN
109: CALL XERBLA( 'ZLA_SYRCOND_X', -INFO )
110: RETURN
111: END IF
112: UP = .FALSE.
113: IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE.
114: *
115: * Compute norm of op(A)*op2(C).
116: *
117: ANORM = 0.0D+0
118: IF ( UP ) THEN
119: DO I = 1, N
120: TMP = 0.0D+0
121: DO J = 1, I
122: TMP = TMP + CABS1( A( J, I ) * X( J ) )
123: END DO
124: DO J = I+1, N
125: TMP = TMP + CABS1( A( I, J ) * X( J ) )
126: END DO
127: RWORK( I ) = TMP
128: ANORM = MAX( ANORM, TMP )
129: END DO
130: ELSE
131: DO I = 1, N
132: TMP = 0.0D+0
133: DO J = 1, I
134: TMP = TMP + CABS1( A( I, J ) * X( J ) )
135: END DO
136: DO J = I+1, N
137: TMP = TMP + CABS1( A( J, I ) * X( J ) )
138: END DO
139: RWORK( I ) = TMP
140: ANORM = MAX( ANORM, TMP )
141: END DO
142: END IF
143: *
144: * Quick return if possible.
145: *
146: IF( N.EQ.0 ) THEN
147: ZLA_SYRCOND_X = 1.0D+0
148: RETURN
149: ELSE IF( ANORM .EQ. 0.0D+0 ) THEN
150: RETURN
151: END IF
152: *
153: * Estimate the norm of inv(op(A)).
154: *
155: AINVNM = 0.0D+0
156: *
157: KASE = 0
158: 10 CONTINUE
159: CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
160: IF( KASE.NE.0 ) THEN
161: IF( KASE.EQ.2 ) THEN
162: *
163: * Multiply by R.
164: *
165: DO I = 1, N
166: WORK( I ) = WORK( I ) * RWORK( I )
167: END DO
168: *
169: IF ( UP ) THEN
170: CALL ZSYTRS( 'U', N, 1, AF, LDAF, IPIV,
171: $ WORK, N, INFO )
172: ELSE
173: CALL ZSYTRS( 'L', N, 1, AF, LDAF, IPIV,
174: $ WORK, N, INFO )
175: ENDIF
176: *
177: * Multiply by inv(X).
178: *
179: DO I = 1, N
180: WORK( I ) = WORK( I ) / X( I )
181: END DO
182: ELSE
183: *
1.5 ! bertrand 184: * Multiply by inv(X**T).
1.1 bertrand 185: *
186: DO I = 1, N
187: WORK( I ) = WORK( I ) / X( I )
188: END DO
189: *
190: IF ( UP ) THEN
191: CALL ZSYTRS( 'U', N, 1, AF, LDAF, IPIV,
192: $ WORK, N, INFO )
193: ELSE
194: CALL ZSYTRS( 'L', N, 1, AF, LDAF, IPIV,
195: $ WORK, N, INFO )
196: END IF
197: *
198: * Multiply by R.
199: *
200: DO I = 1, N
201: WORK( I ) = WORK( I ) * RWORK( I )
202: END DO
203: END IF
204: GO TO 10
205: END IF
206: *
207: * Compute the estimate of the reciprocal condition number.
208: *
209: IF( AINVNM .NE. 0.0D+0 )
210: $ ZLA_SYRCOND_X = 1.0D+0 / AINVNM
211: *
212: RETURN
213: *
214: END