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