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[in] WORK
137: *> \verbatim
138: *> WORK is COMPLEX*16 array, dimension (2*N).
139: *> Workspace.
140: *> \endverbatim
141: *>
142: *> \param[in] 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: *> \date December 2016
157: *
158: *> \ingroup complex16GBcomputational
159: *
160: * =====================================================================
161: DOUBLE PRECISION FUNCTION ZLA_GBRCOND_C( TRANS, N, KL, KU, AB,
162: $ LDAB, AFB, LDAFB, IPIV,
163: $ C, CAPPLY, INFO, WORK,
164: $ RWORK )
165: *
166: * -- LAPACK computational routine (version 3.7.0) --
167: * -- LAPACK is a software package provided by Univ. of Tennessee, --
168: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
169: * December 2016
170: *
171: * .. Scalar Arguments ..
172: CHARACTER TRANS
173: LOGICAL CAPPLY
174: INTEGER N, KL, KU, KD, KE, LDAB, LDAFB, INFO
175: * ..
176: * .. Array Arguments ..
177: INTEGER IPIV( * )
178: COMPLEX*16 AB( LDAB, * ), AFB( LDAFB, * ), WORK( * )
179: DOUBLE PRECISION C( * ), RWORK( * )
180: *
181: *
182: * =====================================================================
183: *
184: * .. Local Scalars ..
185: LOGICAL NOTRANS
186: INTEGER KASE, I, J
187: DOUBLE PRECISION AINVNM, ANORM, TMP
188: COMPLEX*16 ZDUM
189: * ..
190: * .. Local Arrays ..
191: INTEGER ISAVE( 3 )
192: * ..
193: * .. External Functions ..
194: LOGICAL LSAME
195: EXTERNAL LSAME
196: * ..
197: * .. External Subroutines ..
198: EXTERNAL ZLACN2, ZGBTRS, XERBLA
199: * ..
200: * .. Intrinsic Functions ..
201: INTRINSIC ABS, MAX
202: * ..
203: * .. Statement Functions ..
204: DOUBLE PRECISION CABS1
205: * ..
206: * .. Statement Function Definitions ..
207: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
208: * ..
209: * .. Executable Statements ..
210: ZLA_GBRCOND_C = 0.0D+0
211: *
212: INFO = 0
213: NOTRANS = LSAME( TRANS, 'N' )
214: IF ( .NOT. NOTRANS .AND. .NOT. LSAME( TRANS, 'T' ) .AND. .NOT.
215: $ LSAME( TRANS, 'C' ) ) THEN
216: INFO = -1
217: ELSE IF( N.LT.0 ) THEN
218: INFO = -2
219: ELSE IF( KL.LT.0 .OR. KL.GT.N-1 ) THEN
220: INFO = -3
221: ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN
222: INFO = -4
223: ELSE IF( LDAB.LT.KL+KU+1 ) THEN
224: INFO = -6
225: ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN
226: INFO = -8
227: END IF
228: IF( INFO.NE.0 ) THEN
229: CALL XERBLA( 'ZLA_GBRCOND_C', -INFO )
230: RETURN
231: END IF
232: *
233: * Compute norm of op(A)*op2(C).
234: *
235: ANORM = 0.0D+0
236: KD = KU + 1
237: KE = KL + 1
238: IF ( NOTRANS ) THEN
239: DO I = 1, N
240: TMP = 0.0D+0
241: IF ( CAPPLY ) THEN
242: DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
243: TMP = TMP + CABS1( AB( KD+I-J, J ) ) / C( J )
244: END DO
245: ELSE
246: DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
247: TMP = TMP + CABS1( AB( KD+I-J, J ) )
248: END DO
249: END IF
250: RWORK( I ) = TMP
251: ANORM = MAX( ANORM, TMP )
252: END DO
253: ELSE
254: DO I = 1, N
255: TMP = 0.0D+0
256: IF ( CAPPLY ) THEN
257: DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
258: TMP = TMP + CABS1( AB( KE-I+J, I ) ) / C( J )
259: END DO
260: ELSE
261: DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
262: TMP = TMP + CABS1( AB( KE-I+J, I ) )
263: END DO
264: END IF
265: RWORK( I ) = TMP
266: ANORM = MAX( ANORM, TMP )
267: END DO
268: END IF
269: *
270: * Quick return if possible.
271: *
272: IF( N.EQ.0 ) THEN
273: ZLA_GBRCOND_C = 1.0D+0
274: RETURN
275: ELSE IF( ANORM .EQ. 0.0D+0 ) THEN
276: RETURN
277: END IF
278: *
279: * Estimate the norm of inv(op(A)).
280: *
281: AINVNM = 0.0D+0
282: *
283: KASE = 0
284: 10 CONTINUE
285: CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
286: IF( KASE.NE.0 ) THEN
287: IF( KASE.EQ.2 ) THEN
288: *
289: * Multiply by R.
290: *
291: DO I = 1, N
292: WORK( I ) = WORK( I ) * RWORK( I )
293: END DO
294: *
295: IF ( NOTRANS ) THEN
296: CALL ZGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB,
297: $ IPIV, WORK, N, INFO )
298: ELSE
299: CALL ZGBTRS( 'Conjugate transpose', N, KL, KU, 1, AFB,
300: $ LDAFB, IPIV, WORK, N, INFO )
301: ENDIF
302: *
303: * Multiply by inv(C).
304: *
305: IF ( CAPPLY ) THEN
306: DO I = 1, N
307: WORK( I ) = WORK( I ) * C( I )
308: END DO
309: END IF
310: ELSE
311: *
312: * Multiply by inv(C**H).
313: *
314: IF ( CAPPLY ) THEN
315: DO I = 1, N
316: WORK( I ) = WORK( I ) * C( I )
317: END DO
318: END IF
319: *
320: IF ( NOTRANS ) THEN
321: CALL ZGBTRS( 'Conjugate transpose', N, KL, KU, 1, AFB,
322: $ LDAFB, IPIV, WORK, N, INFO )
323: ELSE
324: CALL ZGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB,
325: $ IPIV, WORK, N, INFO )
326: END IF
327: *
328: * Multiply by R.
329: *
330: DO I = 1, N
331: WORK( I ) = WORK( I ) * RWORK( I )
332: END DO
333: END IF
334: GO TO 10
335: END IF
336: *
337: * Compute the estimate of the reciprocal condition number.
338: *
339: IF( AINVNM .NE. 0.0D+0 )
340: $ ZLA_GBRCOND_C = 1.0D+0 / AINVNM
341: *
342: RETURN
343: *
344: 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