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