1: SUBROUTINE DGBRFS( TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB,
2: $ IPIV, B, LDB, X, LDX, FERR, BERR, WORK, IWORK,
3: $ INFO )
4: *
5: * -- LAPACK routine (version 3.2) --
6: * -- LAPACK is a software package provided by Univ. of Tennessee, --
7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8: * November 2006
9: *
10: * Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH.
11: *
12: * .. Scalar Arguments ..
13: CHARACTER TRANS
14: INTEGER INFO, KL, KU, LDAB, LDAFB, LDB, LDX, N, NRHS
15: * ..
16: * .. Array Arguments ..
17: INTEGER IPIV( * ), IWORK( * )
18: DOUBLE PRECISION AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
19: $ BERR( * ), FERR( * ), WORK( * ), X( LDX, * )
20: * ..
21: *
22: * Purpose
23: * =======
24: *
25: * DGBRFS improves the computed solution to a system of linear
26: * equations when the coefficient matrix is banded, and provides
27: * error bounds and backward error estimates for the solution.
28: *
29: * Arguments
30: * =========
31: *
32: * TRANS (input) CHARACTER*1
33: * Specifies the form of the system of equations:
34: * = 'N': A * X = B (No transpose)
35: * = 'T': A**T * X = B (Transpose)
36: * = 'C': A**H * X = B (Conjugate transpose = Transpose)
37: *
38: * N (input) INTEGER
39: * The order of the matrix A. N >= 0.
40: *
41: * KL (input) INTEGER
42: * The number of subdiagonals within the band of A. KL >= 0.
43: *
44: * KU (input) INTEGER
45: * The number of superdiagonals within the band of A. KU >= 0.
46: *
47: * NRHS (input) INTEGER
48: * The number of right hand sides, i.e., the number of columns
49: * of the matrices B and X. NRHS >= 0.
50: *
51: * AB (input) DOUBLE PRECISION array, dimension (LDAB,N)
52: * The original band matrix A, stored 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) DOUBLE PRECISION array, dimension (LDAFB,N)
61: * Details of the LU factorization of the band matrix A, as
62: * computed by DGBTRF. U is stored as an upper triangular band
63: * matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
64: * the multipliers used during the factorization are stored in
65: * 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 DGBTRF; for 1<=i<=N, row i of the
72: * matrix was interchanged with row IPIV(i).
73: *
74: * B (input) DOUBLE PRECISION array, dimension (LDB,NRHS)
75: * The right hand side matrix B.
76: *
77: * LDB (input) INTEGER
78: * The leading dimension of the array B. LDB >= max(1,N).
79: *
80: * X (input/output) DOUBLE PRECISION array, dimension (LDX,NRHS)
81: * On entry, the solution matrix X, as computed by DGBTRS.
82: * On exit, the improved solution matrix X.
83: *
84: * LDX (input) INTEGER
85: * The leading dimension of the array X. LDX >= max(1,N).
86: *
87: * FERR (output) DOUBLE PRECISION array, dimension (NRHS)
88: * The estimated forward error bound for each solution vector
89: * X(j) (the j-th column of the solution matrix X).
90: * If XTRUE is the true solution corresponding to X(j), FERR(j)
91: * is an estimated upper bound for the magnitude of the largest
92: * element in (X(j) - XTRUE) divided by the magnitude of the
93: * largest element in X(j). The estimate is as reliable as
94: * the estimate for RCOND, and is almost always a slight
95: * overestimate of the true error.
96: *
97: * BERR (output) DOUBLE PRECISION array, dimension (NRHS)
98: * The componentwise relative backward error of each solution
99: * vector X(j) (i.e., the smallest relative change in
100: * any element of A or B that makes X(j) an exact solution).
101: *
102: * WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
103: *
104: * IWORK (workspace) INTEGER array, dimension (N)
105: *
106: * INFO (output) INTEGER
107: * = 0: successful exit
108: * < 0: if INFO = -i, the i-th argument had an illegal value
109: *
110: * Internal Parameters
111: * ===================
112: *
113: * ITMAX is the maximum number of steps of iterative refinement.
114: *
115: * =====================================================================
116: *
117: * .. Parameters ..
118: INTEGER ITMAX
119: PARAMETER ( ITMAX = 5 )
120: DOUBLE PRECISION ZERO
121: PARAMETER ( ZERO = 0.0D+0 )
122: DOUBLE PRECISION ONE
123: PARAMETER ( ONE = 1.0D+0 )
124: DOUBLE PRECISION TWO
125: PARAMETER ( TWO = 2.0D+0 )
126: DOUBLE PRECISION THREE
127: PARAMETER ( THREE = 3.0D+0 )
128: * ..
129: * .. Local Scalars ..
130: LOGICAL NOTRAN
131: CHARACTER TRANST
132: INTEGER COUNT, I, J, K, KASE, KK, NZ
133: DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
134: * ..
135: * .. Local Arrays ..
136: INTEGER ISAVE( 3 )
137: * ..
138: * .. External Subroutines ..
139: EXTERNAL DAXPY, DCOPY, DGBMV, DGBTRS, DLACN2, XERBLA
140: * ..
141: * .. Intrinsic Functions ..
142: INTRINSIC ABS, MAX, MIN
143: * ..
144: * .. External Functions ..
145: LOGICAL LSAME
146: DOUBLE PRECISION DLAMCH
147: EXTERNAL LSAME, DLAMCH
148: * ..
149: * .. Executable Statements ..
150: *
151: * Test the input parameters.
152: *
153: INFO = 0
154: NOTRAN = LSAME( TRANS, 'N' )
155: IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
156: $ LSAME( TRANS, 'C' ) ) THEN
157: INFO = -1
158: ELSE IF( N.LT.0 ) THEN
159: INFO = -2
160: ELSE IF( KL.LT.0 ) THEN
161: INFO = -3
162: ELSE IF( KU.LT.0 ) THEN
163: INFO = -4
164: ELSE IF( NRHS.LT.0 ) THEN
165: INFO = -5
166: ELSE IF( LDAB.LT.KL+KU+1 ) THEN
167: INFO = -7
168: ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN
169: INFO = -9
170: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
171: INFO = -12
172: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
173: INFO = -14
174: END IF
175: IF( INFO.NE.0 ) THEN
176: CALL XERBLA( 'DGBRFS', -INFO )
177: RETURN
178: END IF
179: *
180: * Quick return if possible
181: *
182: IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
183: DO 10 J = 1, NRHS
184: FERR( J ) = ZERO
185: BERR( J ) = ZERO
186: 10 CONTINUE
187: RETURN
188: END IF
189: *
190: IF( NOTRAN ) THEN
191: TRANST = 'T'
192: ELSE
193: TRANST = 'N'
194: END IF
195: *
196: * NZ = maximum number of nonzero elements in each row of A, plus 1
197: *
198: NZ = MIN( KL+KU+2, N+1 )
199: EPS = DLAMCH( 'Epsilon' )
200: SAFMIN = DLAMCH( 'Safe minimum' )
201: SAFE1 = NZ*SAFMIN
202: SAFE2 = SAFE1 / EPS
203: *
204: * Do for each right hand side
205: *
206: DO 140 J = 1, NRHS
207: *
208: COUNT = 1
209: LSTRES = THREE
210: 20 CONTINUE
211: *
212: * Loop until stopping criterion is satisfied.
213: *
214: * Compute residual R = B - op(A) * X,
215: * where op(A) = A, A**T, or A**H, depending on TRANS.
216: *
217: CALL DCOPY( N, B( 1, J ), 1, WORK( N+1 ), 1 )
218: CALL DGBMV( TRANS, N, N, KL, KU, -ONE, AB, LDAB, X( 1, J ), 1,
219: $ ONE, WORK( N+1 ), 1 )
220: *
221: * Compute componentwise relative backward error from formula
222: *
223: * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
224: *
225: * where abs(Z) is the componentwise absolute value of the matrix
226: * or vector Z. If the i-th component of the denominator is less
227: * than SAFE2, then SAFE1 is added to the i-th components of the
228: * numerator and denominator before dividing.
229: *
230: DO 30 I = 1, N
231: WORK( I ) = ABS( B( I, J ) )
232: 30 CONTINUE
233: *
234: * Compute abs(op(A))*abs(X) + abs(B).
235: *
236: IF( NOTRAN ) THEN
237: DO 50 K = 1, N
238: KK = KU + 1 - K
239: XK = ABS( X( K, J ) )
240: DO 40 I = MAX( 1, K-KU ), MIN( N, K+KL )
241: WORK( I ) = WORK( I ) + ABS( AB( KK+I, K ) )*XK
242: 40 CONTINUE
243: 50 CONTINUE
244: ELSE
245: DO 70 K = 1, N
246: S = ZERO
247: KK = KU + 1 - K
248: DO 60 I = MAX( 1, K-KU ), MIN( N, K+KL )
249: S = S + ABS( AB( KK+I, K ) )*ABS( X( I, J ) )
250: 60 CONTINUE
251: WORK( K ) = WORK( K ) + S
252: 70 CONTINUE
253: END IF
254: S = ZERO
255: DO 80 I = 1, N
256: IF( WORK( I ).GT.SAFE2 ) THEN
257: S = MAX( S, ABS( WORK( N+I ) ) / WORK( I ) )
258: ELSE
259: S = MAX( S, ( ABS( WORK( N+I ) )+SAFE1 ) /
260: $ ( WORK( I )+SAFE1 ) )
261: END IF
262: 80 CONTINUE
263: BERR( J ) = S
264: *
265: * Test stopping criterion. Continue iterating if
266: * 1) The residual BERR(J) is larger than machine epsilon, and
267: * 2) BERR(J) decreased by at least a factor of 2 during the
268: * last iteration, and
269: * 3) At most ITMAX iterations tried.
270: *
271: IF( BERR( J ).GT.EPS .AND. TWO*BERR( J ).LE.LSTRES .AND.
272: $ COUNT.LE.ITMAX ) THEN
273: *
274: * Update solution and try again.
275: *
276: CALL DGBTRS( TRANS, N, KL, KU, 1, AFB, LDAFB, IPIV,
277: $ WORK( N+1 ), N, INFO )
278: CALL DAXPY( N, ONE, WORK( N+1 ), 1, X( 1, J ), 1 )
279: LSTRES = BERR( J )
280: COUNT = COUNT + 1
281: GO TO 20
282: END IF
283: *
284: * Bound error from formula
285: *
286: * norm(X - XTRUE) / norm(X) .le. FERR =
287: * norm( abs(inv(op(A)))*
288: * ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
289: *
290: * where
291: * norm(Z) is the magnitude of the largest component of Z
292: * inv(op(A)) is the inverse of op(A)
293: * abs(Z) is the componentwise absolute value of the matrix or
294: * vector Z
295: * NZ is the maximum number of nonzeros in any row of A, plus 1
296: * EPS is machine epsilon
297: *
298: * The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
299: * is incremented by SAFE1 if the i-th component of
300: * abs(op(A))*abs(X) + abs(B) is less than SAFE2.
301: *
302: * Use DLACN2 to estimate the infinity-norm of the matrix
303: * inv(op(A)) * diag(W),
304: * where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
305: *
306: DO 90 I = 1, N
307: IF( WORK( I ).GT.SAFE2 ) THEN
308: WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I )
309: ELSE
310: WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I ) + SAFE1
311: END IF
312: 90 CONTINUE
313: *
314: KASE = 0
315: 100 CONTINUE
316: CALL DLACN2( N, WORK( 2*N+1 ), WORK( N+1 ), IWORK, FERR( J ),
317: $ KASE, ISAVE )
318: IF( KASE.NE.0 ) THEN
319: IF( KASE.EQ.1 ) THEN
320: *
321: * Multiply by diag(W)*inv(op(A)**T).
322: *
323: CALL DGBTRS( TRANST, N, KL, KU, 1, AFB, LDAFB, IPIV,
324: $ WORK( N+1 ), N, INFO )
325: DO 110 I = 1, N
326: WORK( N+I ) = WORK( N+I )*WORK( I )
327: 110 CONTINUE
328: ELSE
329: *
330: * Multiply by inv(op(A))*diag(W).
331: *
332: DO 120 I = 1, N
333: WORK( N+I ) = WORK( N+I )*WORK( I )
334: 120 CONTINUE
335: CALL DGBTRS( TRANS, N, KL, KU, 1, AFB, LDAFB, IPIV,
336: $ WORK( N+1 ), N, INFO )
337: END IF
338: GO TO 100
339: END IF
340: *
341: * Normalize error.
342: *
343: LSTRES = ZERO
344: DO 130 I = 1, N
345: LSTRES = MAX( LSTRES, ABS( X( I, J ) ) )
346: 130 CONTINUE
347: IF( LSTRES.NE.ZERO )
348: $ FERR( J ) = FERR( J ) / LSTRES
349: *
350: 140 CONTINUE
351: *
352: RETURN
353: *
354: * End of DGBRFS
355: *
356: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgbrfs.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack