LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dgbrfs.f
Go to the documentation of this file.
1 *> \brief \b DGBRFS
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DGBRFS + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgbrfs.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgbrfs.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgbrfs.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DGBRFS( TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB,
22 * IPIV, B, LDB, X, LDX, FERR, BERR, WORK, IWORK,
23 * INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER TRANS
27 * INTEGER INFO, KL, KU, LDAB, LDAFB, LDB, LDX, N, NRHS
28 * ..
29 * .. Array Arguments ..
30 * INTEGER IPIV( * ), IWORK( * )
31 * DOUBLE PRECISION AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
32 * $ BERR( * ), FERR( * ), WORK( * ), X( LDX, * )
33 * ..
34 *
35 *
36 *> \par Purpose:
37 * =============
38 *>
39 *> \verbatim
40 *>
41 *> DGBRFS improves the computed solution to a system of linear
42 *> equations when the coefficient matrix is banded, and provides
43 *> error bounds and backward error estimates for the solution.
44 *> \endverbatim
45 *
46 * Arguments:
47 * ==========
48 *
49 *> \param[in] TRANS
50 *> \verbatim
51 *> TRANS is CHARACTER*1
52 *> Specifies the form of the system of equations:
53 *> = 'N': A * X = B (No transpose)
54 *> = 'T': A**T * X = B (Transpose)
55 *> = 'C': A**H * X = B (Conjugate transpose = Transpose)
56 *> \endverbatim
57 *>
58 *> \param[in] N
59 *> \verbatim
60 *> N is INTEGER
61 *> The order of the matrix A. N >= 0.
62 *> \endverbatim
63 *>
64 *> \param[in] KL
65 *> \verbatim
66 *> KL is INTEGER
67 *> The number of subdiagonals within the band of A. KL >= 0.
68 *> \endverbatim
69 *>
70 *> \param[in] KU
71 *> \verbatim
72 *> KU is INTEGER
73 *> The number of superdiagonals within the band of A. KU >= 0.
74 *> \endverbatim
75 *>
76 *> \param[in] NRHS
77 *> \verbatim
78 *> NRHS is INTEGER
79 *> The number of right hand sides, i.e., the number of columns
80 *> of the matrices B and X. NRHS >= 0.
81 *> \endverbatim
82 *>
83 *> \param[in] AB
84 *> \verbatim
85 *> AB is DOUBLE PRECISION array, dimension (LDAB,N)
86 *> The original band matrix A, stored in rows 1 to KL+KU+1.
87 *> The j-th column of A is stored in the j-th column of the
88 *> array AB as follows:
89 *> AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(n,j+kl).
90 *> \endverbatim
91 *>
92 *> \param[in] LDAB
93 *> \verbatim
94 *> LDAB is INTEGER
95 *> The leading dimension of the array AB. LDAB >= KL+KU+1.
96 *> \endverbatim
97 *>
98 *> \param[in] AFB
99 *> \verbatim
100 *> AFB is DOUBLE PRECISION array, dimension (LDAFB,N)
101 *> Details of the LU factorization of the band matrix A, as
102 *> computed by DGBTRF. U is stored as an upper triangular band
103 *> matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
104 *> the multipliers used during the factorization are stored in
105 *> rows KL+KU+2 to 2*KL+KU+1.
106 *> \endverbatim
107 *>
108 *> \param[in] LDAFB
109 *> \verbatim
110 *> LDAFB is INTEGER
111 *> The leading dimension of the array AFB. LDAFB >= 2*KL*KU+1.
112 *> \endverbatim
113 *>
114 *> \param[in] IPIV
115 *> \verbatim
116 *> IPIV is INTEGER array, dimension (N)
117 *> The pivot indices from DGBTRF; for 1<=i<=N, row i of the
118 *> matrix was interchanged with row IPIV(i).
119 *> \endverbatim
120 *>
121 *> \param[in] B
122 *> \verbatim
123 *> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
124 *> The right hand side matrix B.
125 *> \endverbatim
126 *>
127 *> \param[in] LDB
128 *> \verbatim
129 *> LDB is INTEGER
130 *> The leading dimension of the array B. LDB >= max(1,N).
131 *> \endverbatim
132 *>
133 *> \param[in,out] X
134 *> \verbatim
135 *> X is DOUBLE PRECISION array, dimension (LDX,NRHS)
136 *> On entry, the solution matrix X, as computed by DGBTRS.
137 *> On exit, the improved solution matrix X.
138 *> \endverbatim
139 *>
140 *> \param[in] LDX
141 *> \verbatim
142 *> LDX is INTEGER
143 *> The leading dimension of the array X. LDX >= max(1,N).
144 *> \endverbatim
145 *>
146 *> \param[out] FERR
147 *> \verbatim
148 *> FERR is DOUBLE PRECISION array, dimension (NRHS)
149 *> The estimated forward error bound for each solution vector
150 *> X(j) (the j-th column of the solution matrix X).
151 *> If XTRUE is the true solution corresponding to X(j), FERR(j)
152 *> is an estimated upper bound for the magnitude of the largest
153 *> element in (X(j) - XTRUE) divided by the magnitude of the
154 *> largest element in X(j). The estimate is as reliable as
155 *> the estimate for RCOND, and is almost always a slight
156 *> overestimate of the true error.
157 *> \endverbatim
158 *>
159 *> \param[out] BERR
160 *> \verbatim
161 *> BERR is DOUBLE PRECISION array, dimension (NRHS)
162 *> The componentwise relative backward error of each solution
163 *> vector X(j) (i.e., the smallest relative change in
164 *> any element of A or B that makes X(j) an exact solution).
165 *> \endverbatim
166 *>
167 *> \param[out] WORK
168 *> \verbatim
169 *> WORK is DOUBLE PRECISION array, dimension (3*N)
170 *> \endverbatim
171 *>
172 *> \param[out] IWORK
173 *> \verbatim
174 *> IWORK is INTEGER array, dimension (N)
175 *> \endverbatim
176 *>
177 *> \param[out] INFO
178 *> \verbatim
179 *> INFO is INTEGER
180 *> = 0: successful exit
181 *> < 0: if INFO = -i, the i-th argument had an illegal value
182 *> \endverbatim
183 *
184 *> \par Internal Parameters:
185 * =========================
186 *>
187 *> \verbatim
188 *> ITMAX is the maximum number of steps of iterative refinement.
189 *> \endverbatim
190 *
191 * Authors:
192 * ========
193 *
194 *> \author Univ. of Tennessee
195 *> \author Univ. of California Berkeley
196 *> \author Univ. of Colorado Denver
197 *> \author NAG Ltd.
198 *
199 *> \date November 2011
200 *
201 *> \ingroup doubleGBcomputational
202 *
203 * =====================================================================
204  SUBROUTINE dgbrfs( TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB,
205  $ ipiv, b, ldb, x, ldx, ferr, berr, work, iwork,
206  $ info )
207 *
208 * -- LAPACK computational routine (version 3.4.0) --
209 * -- LAPACK is a software package provided by Univ. of Tennessee, --
210 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
211 * November 2011
212 *
213 * .. Scalar Arguments ..
214  CHARACTER trans
215  INTEGER info, kl, ku, ldab, ldafb, ldb, ldx, n, nrhs
216 * ..
217 * .. Array Arguments ..
218  INTEGER ipiv( * ), iwork( * )
219  DOUBLE PRECISION ab( ldab, * ), afb( ldafb, * ), b( ldb, * ),
220  $ berr( * ), ferr( * ), work( * ), x( ldx, * )
221 * ..
222 *
223 * =====================================================================
224 *
225 * .. Parameters ..
226  INTEGER itmax
227  parameter( itmax = 5 )
228  DOUBLE PRECISION zero
229  parameter( zero = 0.0d+0 )
230  DOUBLE PRECISION one
231  parameter( one = 1.0d+0 )
232  DOUBLE PRECISION two
233  parameter( two = 2.0d+0 )
234  DOUBLE PRECISION three
235  parameter( three = 3.0d+0 )
236 * ..
237 * .. Local Scalars ..
238  LOGICAL notran
239  CHARACTER transt
240  INTEGER count, i, j, k, kase, kk, nz
241  DOUBLE PRECISION eps, lstres, s, safe1, safe2, safmin, xk
242 * ..
243 * .. Local Arrays ..
244  INTEGER isave( 3 )
245 * ..
246 * .. External Subroutines ..
247  EXTERNAL daxpy, dcopy, dgbmv, dgbtrs, dlacn2, xerbla
248 * ..
249 * .. Intrinsic Functions ..
250  INTRINSIC abs, max, min
251 * ..
252 * .. External Functions ..
253  LOGICAL lsame
254  DOUBLE PRECISION dlamch
255  EXTERNAL lsame, dlamch
256 * ..
257 * .. Executable Statements ..
258 *
259 * Test the input parameters.
260 *
261  info = 0
262  notran = lsame( trans, 'N' )
263  IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
264  $ lsame( trans, 'C' ) ) THEN
265  info = -1
266  ELSE IF( n.LT.0 ) THEN
267  info = -2
268  ELSE IF( kl.LT.0 ) THEN
269  info = -3
270  ELSE IF( ku.LT.0 ) THEN
271  info = -4
272  ELSE IF( nrhs.LT.0 ) THEN
273  info = -5
274  ELSE IF( ldab.LT.kl+ku+1 ) THEN
275  info = -7
276  ELSE IF( ldafb.LT.2*kl+ku+1 ) THEN
277  info = -9
278  ELSE IF( ldb.LT.max( 1, n ) ) THEN
279  info = -12
280  ELSE IF( ldx.LT.max( 1, n ) ) THEN
281  info = -14
282  END IF
283  IF( info.NE.0 ) THEN
284  CALL xerbla( 'DGBRFS', -info )
285  return
286  END IF
287 *
288 * Quick return if possible
289 *
290  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
291  DO 10 j = 1, nrhs
292  ferr( j ) = zero
293  berr( j ) = zero
294  10 continue
295  return
296  END IF
297 *
298  IF( notran ) THEN
299  transt = 'T'
300  ELSE
301  transt = 'N'
302  END IF
303 *
304 * NZ = maximum number of nonzero elements in each row of A, plus 1
305 *
306  nz = min( kl+ku+2, n+1 )
307  eps = dlamch( 'Epsilon' )
308  safmin = dlamch( 'Safe minimum' )
309  safe1 = nz*safmin
310  safe2 = safe1 / eps
311 *
312 * Do for each right hand side
313 *
314  DO 140 j = 1, nrhs
315 *
316  count = 1
317  lstres = three
318  20 continue
319 *
320 * Loop until stopping criterion is satisfied.
321 *
322 * Compute residual R = B - op(A) * X,
323 * where op(A) = A, A**T, or A**H, depending on TRANS.
324 *
325  CALL dcopy( n, b( 1, j ), 1, work( n+1 ), 1 )
326  CALL dgbmv( trans, n, n, kl, ku, -one, ab, ldab, x( 1, j ), 1,
327  $ one, work( n+1 ), 1 )
328 *
329 * Compute componentwise relative backward error from formula
330 *
331 * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
332 *
333 * where abs(Z) is the componentwise absolute value of the matrix
334 * or vector Z. If the i-th component of the denominator is less
335 * than SAFE2, then SAFE1 is added to the i-th components of the
336 * numerator and denominator before dividing.
337 *
338  DO 30 i = 1, n
339  work( i ) = abs( b( i, j ) )
340  30 continue
341 *
342 * Compute abs(op(A))*abs(X) + abs(B).
343 *
344  IF( notran ) THEN
345  DO 50 k = 1, n
346  kk = ku + 1 - k
347  xk = abs( x( k, j ) )
348  DO 40 i = max( 1, k-ku ), min( n, k+kl )
349  work( i ) = work( i ) + abs( ab( kk+i, k ) )*xk
350  40 continue
351  50 continue
352  ELSE
353  DO 70 k = 1, n
354  s = zero
355  kk = ku + 1 - k
356  DO 60 i = max( 1, k-ku ), min( n, k+kl )
357  s = s + abs( ab( kk+i, k ) )*abs( x( i, j ) )
358  60 continue
359  work( k ) = work( k ) + s
360  70 continue
361  END IF
362  s = zero
363  DO 80 i = 1, n
364  IF( work( i ).GT.safe2 ) THEN
365  s = max( s, abs( work( n+i ) ) / work( i ) )
366  ELSE
367  s = max( s, ( abs( work( n+i ) )+safe1 ) /
368  $ ( work( i )+safe1 ) )
369  END IF
370  80 continue
371  berr( j ) = s
372 *
373 * Test stopping criterion. Continue iterating if
374 * 1) The residual BERR(J) is larger than machine epsilon, and
375 * 2) BERR(J) decreased by at least a factor of 2 during the
376 * last iteration, and
377 * 3) At most ITMAX iterations tried.
378 *
379  IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
380  $ count.LE.itmax ) THEN
381 *
382 * Update solution and try again.
383 *
384  CALL dgbtrs( trans, n, kl, ku, 1, afb, ldafb, ipiv,
385  $ work( n+1 ), n, info )
386  CALL daxpy( n, one, work( n+1 ), 1, x( 1, j ), 1 )
387  lstres = berr( j )
388  count = count + 1
389  go to 20
390  END IF
391 *
392 * Bound error from formula
393 *
394 * norm(X - XTRUE) / norm(X) .le. FERR =
395 * norm( abs(inv(op(A)))*
396 * ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
397 *
398 * where
399 * norm(Z) is the magnitude of the largest component of Z
400 * inv(op(A)) is the inverse of op(A)
401 * abs(Z) is the componentwise absolute value of the matrix or
402 * vector Z
403 * NZ is the maximum number of nonzeros in any row of A, plus 1
404 * EPS is machine epsilon
405 *
406 * The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
407 * is incremented by SAFE1 if the i-th component of
408 * abs(op(A))*abs(X) + abs(B) is less than SAFE2.
409 *
410 * Use DLACN2 to estimate the infinity-norm of the matrix
411 * inv(op(A)) * diag(W),
412 * where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
413 *
414  DO 90 i = 1, n
415  IF( work( i ).GT.safe2 ) THEN
416  work( i ) = abs( work( n+i ) ) + nz*eps*work( i )
417  ELSE
418  work( i ) = abs( work( n+i ) ) + nz*eps*work( i ) + safe1
419  END IF
420  90 continue
421 *
422  kase = 0
423  100 continue
424  CALL dlacn2( n, work( 2*n+1 ), work( n+1 ), iwork, ferr( j ),
425  $ kase, isave )
426  IF( kase.NE.0 ) THEN
427  IF( kase.EQ.1 ) THEN
428 *
429 * Multiply by diag(W)*inv(op(A)**T).
430 *
431  CALL dgbtrs( transt, n, kl, ku, 1, afb, ldafb, ipiv,
432  $ work( n+1 ), n, info )
433  DO 110 i = 1, n
434  work( n+i ) = work( n+i )*work( i )
435  110 continue
436  ELSE
437 *
438 * Multiply by inv(op(A))*diag(W).
439 *
440  DO 120 i = 1, n
441  work( n+i ) = work( n+i )*work( i )
442  120 continue
443  CALL dgbtrs( trans, n, kl, ku, 1, afb, ldafb, ipiv,
444  $ work( n+1 ), n, info )
445  END IF
446  go to 100
447  END IF
448 *
449 * Normalize error.
450 *
451  lstres = zero
452  DO 130 i = 1, n
453  lstres = max( lstres, abs( x( i, j ) ) )
454  130 continue
455  IF( lstres.NE.zero )
456  $ ferr( j ) = ferr( j ) / lstres
457 *
458  140 continue
459 *
460  return
461 *
462 * End of DGBRFS
463 *
464  END