![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zla_gercond_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_gercond_x.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: DOUBLE PRECISION FUNCTION ZLA_GERCOND_X( TRANS, 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 TRANS
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_GERCOND_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: * TRANS (input) CHARACTER*1
35: * Specifies the form of the system of equations:
36: * = 'N': A * X = B (No transpose)
37: * = 'T': A**T * X = B (Transpose)
38: * = 'C': A**H * X = B (Conjugate Transpose = Transpose)
39: *
40: * N (input) INTEGER
41: * The number of linear equations, i.e., the order of the
42: * matrix A. N >= 0.
43: *
44: * A (input) COMPLEX*16 array, dimension (LDA,N)
45: * On entry, the N-by-N matrix A.
46: *
47: * LDA (input) INTEGER
48: * The leading dimension of the array A. LDA >= max(1,N).
49: *
50: * AF (input) COMPLEX*16 array, dimension (LDAF,N)
51: * The factors L and U from the factorization
52: * A = P*L*U as computed by ZGETRF.
53: *
54: * LDAF (input) INTEGER
55: * The leading dimension of the array AF. LDAF >= max(1,N).
56: *
57: * IPIV (input) INTEGER array, dimension (N)
58: * The pivot indices from the factorization A = P*L*U
59: * as computed by ZGETRF; row i of the matrix was interchanged
60: * with row IPIV(i).
61: *
62: * X (input) COMPLEX*16 array, dimension (N)
63: * The vector X in the formula op(A) * diag(X).
64: *
65: * INFO (output) INTEGER
66: * = 0: Successful exit.
67: * i > 0: The ith argument is invalid.
68: *
69: * WORK (input) COMPLEX*16 array, dimension (2*N).
70: * Workspace.
71: *
72: * RWORK (input) DOUBLE PRECISION array, dimension (N).
73: * Workspace.
74: *
75: * =====================================================================
76: *
77: * .. Local Scalars ..
78: LOGICAL NOTRANS
79: INTEGER KASE
80: DOUBLE PRECISION AINVNM, ANORM, TMP
81: INTEGER I, J
82: COMPLEX*16 ZDUM
83: * ..
84: * .. Local Arrays ..
85: INTEGER ISAVE( 3 )
86: * ..
87: * .. External Functions ..
88: LOGICAL LSAME
89: EXTERNAL LSAME
90: * ..
91: * .. External Subroutines ..
92: EXTERNAL ZLACN2, ZGETRS, XERBLA
93: * ..
94: * .. Intrinsic Functions ..
95: INTRINSIC ABS, MAX, REAL, DIMAG
96: * ..
97: * .. Statement Functions ..
98: DOUBLE PRECISION CABS1
99: * ..
100: * .. Statement Function Definitions ..
101: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
102: * ..
103: * .. Executable Statements ..
104: *
105: ZLA_GERCOND_X = 0.0D+0
106: *
107: INFO = 0
108: NOTRANS = LSAME( TRANS, 'N' )
109: IF ( .NOT. NOTRANS .AND. .NOT. LSAME( TRANS, 'T' ) .AND. .NOT.
110: $ LSAME( TRANS, 'C' ) ) THEN
111: INFO = -1
112: ELSE IF( N.LT.0 ) THEN
113: INFO = -2
114: END IF
115: IF( INFO.NE.0 ) THEN
116: CALL XERBLA( 'ZLA_GERCOND_X', -INFO )
117: RETURN
118: END IF
119: *
120: * Compute norm of op(A)*op2(C).
121: *
122: ANORM = 0.0D+0
123: IF ( NOTRANS ) THEN
124: DO I = 1, N
125: TMP = 0.0D+0
126: DO J = 1, N
127: TMP = TMP + CABS1( A( I, J ) * X( J ) )
128: END DO
129: RWORK( I ) = TMP
130: ANORM = MAX( ANORM, TMP )
131: END DO
132: ELSE
133: DO I = 1, N
134: TMP = 0.0D+0
135: DO J = 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_GERCOND_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: * Multiply by R.
162: DO I = 1, N
163: WORK( I ) = WORK( I ) * RWORK( I )
164: END DO
165: *
166: IF ( NOTRANS ) THEN
167: CALL ZGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
168: $ WORK, N, INFO )
169: ELSE
170: CALL ZGETRS( 'Conjugate transpose', N, 1, AF, LDAF, IPIV,
171: $ WORK, N, INFO )
172: ENDIF
173: *
174: * Multiply by inv(X).
175: *
176: DO I = 1, N
177: WORK( I ) = WORK( I ) / X( I )
178: END DO
179: ELSE
180: *
181: * Multiply by inv(X').
182: *
183: DO I = 1, N
184: WORK( I ) = WORK( I ) / X( I )
185: END DO
186: *
187: IF ( NOTRANS ) THEN
188: CALL ZGETRS( 'Conjugate transpose', N, 1, AF, LDAF, IPIV,
189: $ WORK, N, INFO )
190: ELSE
191: CALL ZGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
192: $ WORK, N, INFO )
193: END IF
194: *
195: * Multiply by R.
196: *
197: DO I = 1, N
198: WORK( I ) = WORK( I ) * RWORK( I )
199: END DO
200: END IF
201: GO TO 10
202: END IF
203: *
204: * Compute the estimate of the reciprocal condition number.
205: *
206: IF( AINVNM .NE. 0.0D+0 )
207: $ ZLA_GERCOND_X = 1.0D+0 / AINVNM
208: *
209: RETURN
210: *
211: END