![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zla_gbrcond_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_gbrcond_c.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: DOUBLE PRECISION FUNCTION ZLA_GBRCOND_C( TRANS, N, KL, KU, AB,
2: $ LDAB, AFB, LDAFB, IPIV,
3: $ C, CAPPLY, INFO, WORK,
4: $ RWORK )
5: *
6: * -- LAPACK routine (version 3.2.1) --
7: * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
8: * -- Jason Riedy of Univ. of California Berkeley. --
9: * -- April 2009 --
10: *
11: * -- LAPACK is a software package provided by Univ. of Tennessee, --
12: * -- Univ. of California Berkeley and NAG Ltd. --
13: *
14: IMPLICIT NONE
15: * ..
16: * .. Scalar Arguments ..
17: CHARACTER TRANS
18: LOGICAL CAPPLY
19: INTEGER N, KL, KU, KD, KE, LDAB, LDAFB, INFO
20: * ..
21: * .. Array Arguments ..
22: INTEGER IPIV( * )
23: COMPLEX*16 AB( LDAB, * ), AFB( LDAFB, * ), WORK( * )
24: DOUBLE PRECISION C( * ), RWORK( * )
25: *
26: *
27: * Purpose
28: * =======
29: *
30: * ZLA_GBRCOND_C Computes the infinity norm condition number of
31: * op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector.
32: *
33: * Arguments
34: * =========
35: *
36: * TRANS (input) CHARACTER*1
37: * Specifies the form of the system of equations:
38: * = 'N': A * X = B (No transpose)
39: * = 'T': A**T * X = B (Transpose)
40: * = 'C': A**H * X = B (Conjugate Transpose = Transpose)
41: *
42: * N (input) INTEGER
43: * The number of linear equations, i.e., the order of the
44: * matrix A. N >= 0.
45: *
46: * KL (input) INTEGER
47: * The number of subdiagonals within the band of A. KL >= 0.
48: *
49: * KU (input) INTEGER
50: * The number of superdiagonals within the band of A. KU >= 0.
51: *
52: * AB (input) COMPLEX*16 array, dimension (LDAB,N)
53: * On entry, the matrix A in band storage, in rows 1 to KL+KU+1.
54: * The j-th column of A is stored in the j-th column of the
55: * array AB as follows:
56: * AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)
57: *
58: * LDAB (input) INTEGER
59: * The leading dimension of the array AB. LDAB >= KL+KU+1.
60: *
61: * AFB (input) COMPLEX*16 array, dimension (LDAFB,N)
62: * Details of the LU factorization of the band matrix A, as
63: * computed by ZGBTRF. U is stored as an upper triangular
64: * band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1,
65: * and the multipliers used during the factorization are stored
66: * in rows KL+KU+2 to 2*KL+KU+1.
67: *
68: * LDAFB (input) INTEGER
69: * The leading dimension of the array AFB. LDAFB >= 2*KL+KU+1.
70: *
71: * IPIV (input) INTEGER array, dimension (N)
72: * The pivot indices from the factorization A = P*L*U
73: * as computed by ZGBTRF; row i of the matrix was interchanged
74: * with row IPIV(i).
75: *
76: * C (input) DOUBLE PRECISION array, dimension (N)
77: * The vector C in the formula op(A) * inv(diag(C)).
78: *
79: * CAPPLY (input) LOGICAL
80: * If .TRUE. then access the vector C in the formula above.
81: *
82: * INFO (output) INTEGER
83: * = 0: Successful exit.
84: * i > 0: The ith argument is invalid.
85: *
86: * WORK (input) COMPLEX*16 array, dimension (2*N).
87: * Workspace.
88: *
89: * RWORK (input) DOUBLE PRECISION array, dimension (N).
90: * Workspace.
91: *
92: * =====================================================================
93: *
94: * .. Local Scalars ..
95: LOGICAL NOTRANS
96: INTEGER KASE, I, J
97: DOUBLE PRECISION AINVNM, ANORM, TMP
98: COMPLEX*16 ZDUM
99: * ..
100: * .. Local Arrays ..
101: INTEGER ISAVE( 3 )
102: * ..
103: * .. External Functions ..
104: LOGICAL LSAME
105: EXTERNAL LSAME
106: * ..
107: * .. External Subroutines ..
108: EXTERNAL ZLACN2, ZGBTRS, XERBLA
109: * ..
110: * .. Intrinsic Functions ..
111: INTRINSIC ABS, MAX
112: * ..
113: * .. Statement Functions ..
114: DOUBLE PRECISION CABS1
115: * ..
116: * .. Statement Function Definitions ..
117: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
118: * ..
119: * .. Executable Statements ..
120: ZLA_GBRCOND_C = 0.0D+0
121: *
122: INFO = 0
123: NOTRANS = LSAME( TRANS, 'N' )
124: IF ( .NOT. NOTRANS .AND. .NOT. LSAME( TRANS, 'T' ) .AND. .NOT.
125: $ LSAME( TRANS, 'C' ) ) THEN
126: INFO = -1
127: ELSE IF( N.LT.0 ) THEN
128: INFO = -2
129: ELSE IF( KL.LT.0 .OR. KL.GT.N-1 ) THEN
130: INFO = -3
131: ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN
132: INFO = -4
133: ELSE IF( LDAB.LT.KL+KU+1 ) THEN
134: INFO = -6
135: ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN
136: INFO = -8
137: END IF
138: IF( INFO.NE.0 ) THEN
139: CALL XERBLA( 'ZLA_GBRCOND_C', -INFO )
140: RETURN
141: END IF
142: *
143: * Compute norm of op(A)*op2(C).
144: *
145: ANORM = 0.0D+0
146: KD = KU + 1
147: KE = KL + 1
148: IF ( NOTRANS ) THEN
149: DO I = 1, N
150: TMP = 0.0D+0
151: IF ( CAPPLY ) THEN
152: DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
153: TMP = TMP + CABS1( AB( KD+I-J, J ) ) / C( J )
154: END DO
155: ELSE
156: DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
157: TMP = TMP + CABS1( AB( KD+I-J, J ) )
158: END DO
159: END IF
160: RWORK( I ) = TMP
161: ANORM = MAX( ANORM, TMP )
162: END DO
163: ELSE
164: DO I = 1, N
165: TMP = 0.0D+0
166: IF ( CAPPLY ) THEN
167: DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
168: TMP = TMP + CABS1( AB( KE-I+J, I ) ) / C( J )
169: END DO
170: ELSE
171: DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
172: TMP = TMP + CABS1( AB( KE-I+J, I ) )
173: END DO
174: END IF
175: RWORK( I ) = TMP
176: ANORM = MAX( ANORM, TMP )
177: END DO
178: END IF
179: *
180: * Quick return if possible.
181: *
182: IF( N.EQ.0 ) THEN
183: ZLA_GBRCOND_C = 1.0D+0
184: RETURN
185: ELSE IF( ANORM .EQ. 0.0D+0 ) THEN
186: RETURN
187: END IF
188: *
189: * Estimate the norm of inv(op(A)).
190: *
191: AINVNM = 0.0D+0
192: *
193: KASE = 0
194: 10 CONTINUE
195: CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
196: IF( KASE.NE.0 ) THEN
197: IF( KASE.EQ.2 ) THEN
198: *
199: * Multiply by R.
200: *
201: DO I = 1, N
202: WORK( I ) = WORK( I ) * RWORK( I )
203: END DO
204: *
205: IF ( NOTRANS ) THEN
206: CALL ZGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB,
207: $ IPIV, WORK, N, INFO )
208: ELSE
209: CALL ZGBTRS( 'Conjugate transpose', N, KL, KU, 1, AFB,
210: $ LDAFB, IPIV, WORK, N, INFO )
211: ENDIF
212: *
213: * Multiply by inv(C).
214: *
215: IF ( CAPPLY ) THEN
216: DO I = 1, N
217: WORK( I ) = WORK( I ) * C( I )
218: END DO
219: END IF
220: ELSE
221: *
222: * Multiply by inv(C').
223: *
224: IF ( CAPPLY ) THEN
225: DO I = 1, N
226: WORK( I ) = WORK( I ) * C( I )
227: END DO
228: END IF
229: *
230: IF ( NOTRANS ) THEN
231: CALL ZGBTRS( 'Conjugate transpose', N, KL, KU, 1, AFB,
232: $ LDAFB, IPIV, WORK, N, INFO )
233: ELSE
234: CALL ZGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB,
235: $ IPIV, WORK, N, INFO )
236: END IF
237: *
238: * Multiply by R.
239: *
240: DO I = 1, N
241: WORK( I ) = WORK( I ) * RWORK( I )
242: END DO
243: END IF
244: GO TO 10
245: END IF
246: *
247: * Compute the estimate of the reciprocal condition number.
248: *
249: IF( AINVNM .NE. 0.0D+0 )
250: $ ZLA_GBRCOND_C = 1.0D+0 / AINVNM
251: *
252: RETURN
253: *
254: END