1: *> \brief \b ZLA_GBRCOND_C computes the infinity norm condition number of op(A)*inv(diag(c)) for general banded matrices.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZLA_GBRCOND_C + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zla_gbrcond_c.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zla_gbrcond_c.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zla_gbrcond_c.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * DOUBLE PRECISION FUNCTION ZLA_GBRCOND_C( TRANS, N, KL, KU, AB,
22: * LDAB, AFB, LDAFB, IPIV,
23: * C, CAPPLY, INFO, WORK,
24: * RWORK )
25: *
26: * .. Scalar Arguments ..
27: * CHARACTER TRANS
28: * LOGICAL CAPPLY
29: * INTEGER N, KL, KU, KD, KE, LDAB, LDAFB, INFO
30: * ..
31: * .. Array Arguments ..
32: * INTEGER IPIV( * )
33: * COMPLEX*16 AB( LDAB, * ), AFB( LDAFB, * ), WORK( * )
34: * DOUBLE PRECISION C( * ), RWORK( * )
35: *
36: *
37: *
38: *> \par Purpose:
39: * =============
40: *>
41: *> \verbatim
42: *>
43: *> ZLA_GBRCOND_C Computes the infinity norm condition number of
44: *> op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector.
45: *> \endverbatim
46: *
47: * Arguments:
48: * ==========
49: *
50: *> \param[in] TRANS
51: *> \verbatim
52: *> TRANS is CHARACTER*1
53: *> Specifies the form of the system of equations:
54: *> = 'N': A * X = B (No transpose)
55: *> = 'T': A**T * X = B (Transpose)
56: *> = 'C': A**H * X = B (Conjugate Transpose = Transpose)
57: *> \endverbatim
58: *>
59: *> \param[in] N
60: *> \verbatim
61: *> N is INTEGER
62: *> The number of linear equations, i.e., the order of the
63: *> matrix A. N >= 0.
64: *> \endverbatim
65: *>
66: *> \param[in] KL
67: *> \verbatim
68: *> KL is INTEGER
69: *> The number of subdiagonals within the band of A. KL >= 0.
70: *> \endverbatim
71: *>
72: *> \param[in] KU
73: *> \verbatim
74: *> KU is INTEGER
75: *> The number of superdiagonals within the band of A. KU >= 0.
76: *> \endverbatim
77: *>
78: *> \param[in] AB
79: *> \verbatim
80: *> AB is COMPLEX*16 array, dimension (LDAB,N)
81: *> On entry, the matrix A in band storage, in rows 1 to KL+KU+1.
82: *> The j-th column of A is stored in the j-th column of the
83: *> array AB as follows:
84: *> AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)
85: *> \endverbatim
86: *>
87: *> \param[in] LDAB
88: *> \verbatim
89: *> LDAB is INTEGER
90: *> The leading dimension of the array AB. LDAB >= KL+KU+1.
91: *> \endverbatim
92: *>
93: *> \param[in] AFB
94: *> \verbatim
95: *> AFB is COMPLEX*16 array, dimension (LDAFB,N)
96: *> Details of the LU factorization of the band matrix A, as
97: *> computed by ZGBTRF. U is stored as an upper triangular
98: *> band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1,
99: *> and the multipliers used during the factorization are stored
100: *> in rows KL+KU+2 to 2*KL+KU+1.
101: *> \endverbatim
102: *>
103: *> \param[in] LDAFB
104: *> \verbatim
105: *> LDAFB is INTEGER
106: *> The leading dimension of the array AFB. LDAFB >= 2*KL+KU+1.
107: *> \endverbatim
108: *>
109: *> \param[in] IPIV
110: *> \verbatim
111: *> IPIV is INTEGER array, dimension (N)
112: *> The pivot indices from the factorization A = P*L*U
113: *> as computed by ZGBTRF; row i of the matrix was interchanged
114: *> with row IPIV(i).
115: *> \endverbatim
116: *>
117: *> \param[in] C
118: *> \verbatim
119: *> C is DOUBLE PRECISION array, dimension (N)
120: *> The vector C in the formula op(A) * inv(diag(C)).
121: *> \endverbatim
122: *>
123: *> \param[in] CAPPLY
124: *> \verbatim
125: *> CAPPLY is LOGICAL
126: *> If .TRUE. then access the vector C in the formula above.
127: *> \endverbatim
128: *>
129: *> \param[out] INFO
130: *> \verbatim
131: *> INFO is INTEGER
132: *> = 0: Successful exit.
133: *> i > 0: The ith argument is invalid.
134: *> \endverbatim
135: *>
136: *> \param[out] WORK
137: *> \verbatim
138: *> WORK is COMPLEX*16 array, dimension (2*N).
139: *> Workspace.
140: *> \endverbatim
141: *>
142: *> \param[out] RWORK
143: *> \verbatim
144: *> RWORK is DOUBLE PRECISION array, dimension (N).
145: *> Workspace.
146: *> \endverbatim
147: *
148: * Authors:
149: * ========
150: *
151: *> \author Univ. of Tennessee
152: *> \author Univ. of California Berkeley
153: *> \author Univ. of Colorado Denver
154: *> \author NAG Ltd.
155: *
156: *> \ingroup complex16GBcomputational
157: *
158: * =====================================================================
159: DOUBLE PRECISION FUNCTION ZLA_GBRCOND_C( TRANS, N, KL, KU, AB,
160: $ LDAB, AFB, LDAFB, IPIV,
161: $ C, CAPPLY, INFO, WORK,
162: $ RWORK )
163: *
164: * -- LAPACK computational routine --
165: * -- LAPACK is a software package provided by Univ. of Tennessee, --
166: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
167: *
168: * .. Scalar Arguments ..
169: CHARACTER TRANS
170: LOGICAL CAPPLY
171: INTEGER N, KL, KU, KD, KE, LDAB, LDAFB, INFO
172: * ..
173: * .. Array Arguments ..
174: INTEGER IPIV( * )
175: COMPLEX*16 AB( LDAB, * ), AFB( LDAFB, * ), WORK( * )
176: DOUBLE PRECISION C( * ), RWORK( * )
177: *
178: *
179: * =====================================================================
180: *
181: * .. Local Scalars ..
182: LOGICAL NOTRANS
183: INTEGER KASE, I, J
184: DOUBLE PRECISION AINVNM, ANORM, TMP
185: COMPLEX*16 ZDUM
186: * ..
187: * .. Local Arrays ..
188: INTEGER ISAVE( 3 )
189: * ..
190: * .. External Functions ..
191: LOGICAL LSAME
192: EXTERNAL LSAME
193: * ..
194: * .. External Subroutines ..
195: EXTERNAL ZLACN2, ZGBTRS, XERBLA
196: * ..
197: * .. Intrinsic Functions ..
198: INTRINSIC ABS, MAX
199: * ..
200: * .. Statement Functions ..
201: DOUBLE PRECISION CABS1
202: * ..
203: * .. Statement Function Definitions ..
204: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
205: * ..
206: * .. Executable Statements ..
207: ZLA_GBRCOND_C = 0.0D+0
208: *
209: INFO = 0
210: NOTRANS = LSAME( TRANS, 'N' )
211: IF ( .NOT. NOTRANS .AND. .NOT. LSAME( TRANS, 'T' ) .AND. .NOT.
212: $ LSAME( TRANS, 'C' ) ) THEN
213: INFO = -1
214: ELSE IF( N.LT.0 ) THEN
215: INFO = -2
216: ELSE IF( KL.LT.0 .OR. KL.GT.N-1 ) THEN
217: INFO = -3
218: ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN
219: INFO = -4
220: ELSE IF( LDAB.LT.KL+KU+1 ) THEN
221: INFO = -6
222: ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN
223: INFO = -8
224: END IF
225: IF( INFO.NE.0 ) THEN
226: CALL XERBLA( 'ZLA_GBRCOND_C', -INFO )
227: RETURN
228: END IF
229: *
230: * Compute norm of op(A)*op2(C).
231: *
232: ANORM = 0.0D+0
233: KD = KU + 1
234: KE = KL + 1
235: IF ( NOTRANS ) THEN
236: DO I = 1, N
237: TMP = 0.0D+0
238: IF ( CAPPLY ) THEN
239: DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
240: TMP = TMP + CABS1( AB( KD+I-J, J ) ) / C( J )
241: END DO
242: ELSE
243: DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
244: TMP = TMP + CABS1( AB( KD+I-J, J ) )
245: END DO
246: END IF
247: RWORK( I ) = TMP
248: ANORM = MAX( ANORM, TMP )
249: END DO
250: ELSE
251: DO I = 1, N
252: TMP = 0.0D+0
253: IF ( CAPPLY ) THEN
254: DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
255: TMP = TMP + CABS1( AB( KE-I+J, I ) ) / C( J )
256: END DO
257: ELSE
258: DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
259: TMP = TMP + CABS1( AB( KE-I+J, I ) )
260: END DO
261: END IF
262: RWORK( I ) = TMP
263: ANORM = MAX( ANORM, TMP )
264: END DO
265: END IF
266: *
267: * Quick return if possible.
268: *
269: IF( N.EQ.0 ) THEN
270: ZLA_GBRCOND_C = 1.0D+0
271: RETURN
272: ELSE IF( ANORM .EQ. 0.0D+0 ) THEN
273: RETURN
274: END IF
275: *
276: * Estimate the norm of inv(op(A)).
277: *
278: AINVNM = 0.0D+0
279: *
280: KASE = 0
281: 10 CONTINUE
282: CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
283: IF( KASE.NE.0 ) THEN
284: IF( KASE.EQ.2 ) THEN
285: *
286: * Multiply by R.
287: *
288: DO I = 1, N
289: WORK( I ) = WORK( I ) * RWORK( I )
290: END DO
291: *
292: IF ( NOTRANS ) THEN
293: CALL ZGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB,
294: $ IPIV, WORK, N, INFO )
295: ELSE
296: CALL ZGBTRS( 'Conjugate transpose', N, KL, KU, 1, AFB,
297: $ LDAFB, IPIV, WORK, N, INFO )
298: ENDIF
299: *
300: * Multiply by inv(C).
301: *
302: IF ( CAPPLY ) THEN
303: DO I = 1, N
304: WORK( I ) = WORK( I ) * C( I )
305: END DO
306: END IF
307: ELSE
308: *
309: * Multiply by inv(C**H).
310: *
311: IF ( CAPPLY ) THEN
312: DO I = 1, N
313: WORK( I ) = WORK( I ) * C( I )
314: END DO
315: END IF
316: *
317: IF ( NOTRANS ) THEN
318: CALL ZGBTRS( 'Conjugate transpose', N, KL, KU, 1, AFB,
319: $ LDAFB, IPIV, WORK, N, INFO )
320: ELSE
321: CALL ZGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB,
322: $ IPIV, WORK, N, INFO )
323: END IF
324: *
325: * Multiply by R.
326: *
327: DO I = 1, N
328: WORK( I ) = WORK( I ) * RWORK( I )
329: END DO
330: END IF
331: GO TO 10
332: END IF
333: *
334: * Compute the estimate of the reciprocal condition number.
335: *
336: IF( AINVNM .NE. 0.0D+0 )
337: $ ZLA_GBRCOND_C = 1.0D+0 / AINVNM
338: *
339: RETURN
340: *
341: * End of ZLA_GBRCOND_C
342: *
343: END
CVSweb interface <joel.bertrand@systella.fr>
![[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