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