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*> \ingroup doubleGBcomputational
200*
201*  =====================================================================
202      SUBROUTINE DGBRFS( TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB,
203     $                   IPIV, B, LDB, X, LDX, FERR, BERR, WORK, IWORK,
204     $                   INFO )
205*
206*  -- LAPACK computational routine --
207*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
208*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
209*
210*     .. Scalar Arguments ..
211      CHARACTER          TRANS
212      INTEGER            INFO, KL, KU, LDAB, LDAFB, LDB, LDX, N, NRHS
213*     ..
214*     .. Array Arguments ..
215      INTEGER            IPIV( * ), IWORK( * )
216      DOUBLE PRECISION   AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
217     $                   BERR( * ), FERR( * ), WORK( * ), X( LDX, * )
218*     ..
219*
220*  =====================================================================
221*
222*     .. Parameters ..
223      INTEGER            ITMAX
224      PARAMETER          ( ITMAX = 5 )
225      DOUBLE PRECISION   ZERO
226      PARAMETER          ( ZERO = 0.0D+0 )
227      DOUBLE PRECISION   ONE
228      PARAMETER          ( ONE = 1.0D+0 )
229      DOUBLE PRECISION   TWO
230      PARAMETER          ( TWO = 2.0D+0 )
231      DOUBLE PRECISION   THREE
232      PARAMETER          ( THREE = 3.0D+0 )
233*     ..
234*     .. Local Scalars ..
235      LOGICAL            NOTRAN
236      CHARACTER          TRANST
237      INTEGER            COUNT, I, J, K, KASE, KK, NZ
238      DOUBLE PRECISION   EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
239*     ..
240*     .. Local Arrays ..
241      INTEGER            ISAVE( 3 )
242*     ..
243*     .. External Subroutines ..
244      EXTERNAL           DAXPY, DCOPY, DGBMV, DGBTRS, DLACN2, XERBLA
245*     ..
246*     .. Intrinsic Functions ..
247      INTRINSIC          ABS, MAX, MIN
248*     ..
249*     .. External Functions ..
250      LOGICAL            LSAME
251      DOUBLE PRECISION   DLAMCH
252      EXTERNAL           LSAME, DLAMCH
253*     ..
254*     .. Executable Statements ..
255*
256*     Test the input parameters.
257*
258      INFO = 0
259      NOTRAN = LSAME( TRANS, 'N' )
260      IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
261     $    LSAME( TRANS, 'C' ) ) THEN
262         INFO = -1
263      ELSE IF( N.LT.0 ) THEN
264         INFO = -2
265      ELSE IF( KL.LT.0 ) THEN
266         INFO = -3
267      ELSE IF( KU.LT.0 ) THEN
268         INFO = -4
269      ELSE IF( NRHS.LT.0 ) THEN
270         INFO = -5
271      ELSE IF( LDAB.LT.KL+KU+1 ) THEN
272         INFO = -7
273      ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN
274         INFO = -9
275      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
276         INFO = -12
277      ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
278         INFO = -14
279      END IF
280      IF( INFO.NE.0 ) THEN
281         CALL XERBLA( 'DGBRFS', -INFO )
282         RETURN
283      END IF
284*
285*     Quick return if possible
286*
287      IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
288         DO 10 J = 1, NRHS
289            FERR( J ) = ZERO
290            BERR( J ) = ZERO
291   10    CONTINUE
292         RETURN
293      END IF
294*
295      IF( NOTRAN ) THEN
296         TRANST = 'T'
297      ELSE
298         TRANST = 'N'
299      END IF
300*
301*     NZ = maximum number of nonzero elements in each row of A, plus 1
302*
303      NZ = MIN( KL+KU+2, N+1 )
304      EPS = DLAMCH( 'Epsilon' )
305      SAFMIN = DLAMCH( 'Safe minimum' )
306      SAFE1 = NZ*SAFMIN
307      SAFE2 = SAFE1 / EPS
308*
309*     Do for each right hand side
310*
311      DO 140 J = 1, NRHS
312*
313         COUNT = 1
314         LSTRES = THREE
315   20    CONTINUE
316*
317*        Loop until stopping criterion is satisfied.
318*
319*        Compute residual R = B - op(A) * X,
320*        where op(A) = A, A**T, or A**H, depending on TRANS.
321*
322         CALL DCOPY( N, B( 1, J ), 1, WORK( N+1 ), 1 )
323         CALL DGBMV( TRANS, N, N, KL, KU, -ONE, AB, LDAB, X( 1, J ), 1,
324     $               ONE, WORK( N+1 ), 1 )
325*
326*        Compute componentwise relative backward error from formula
327*
328*        max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
329*
330*        where abs(Z) is the componentwise absolute value of the matrix
331*        or vector Z.  If the i-th component of the denominator is less
332*        than SAFE2, then SAFE1 is added to the i-th components of the
333*        numerator and denominator before dividing.
334*
335         DO 30 I = 1, N
336            WORK( I ) = ABS( B( I, J ) )
337   30    CONTINUE
338*
339*        Compute abs(op(A))*abs(X) + abs(B).
340*
341         IF( NOTRAN ) THEN
342            DO 50 K = 1, N
343               KK = KU + 1 - K
344               XK = ABS( X( K, J ) )
345               DO 40 I = MAX( 1, K-KU ), MIN( N, K+KL )
346                  WORK( I ) = WORK( I ) + ABS( AB( KK+I, K ) )*XK
347   40          CONTINUE
348   50       CONTINUE
349         ELSE
350            DO 70 K = 1, N
351               S = ZERO
352               KK = KU + 1 - K
353               DO 60 I = MAX( 1, K-KU ), MIN( N, K+KL )
354                  S = S + ABS( AB( KK+I, K ) )*ABS( X( I, J ) )
355   60          CONTINUE
356               WORK( K ) = WORK( K ) + S
357   70       CONTINUE
358         END IF
359         S = ZERO
360         DO 80 I = 1, N
361            IF( WORK( I ).GT.SAFE2 ) THEN
362               S = MAX( S, ABS( WORK( N+I ) ) / WORK( I ) )
363            ELSE
364               S = MAX( S, ( ABS( WORK( N+I ) )+SAFE1 ) /
365     $             ( WORK( I )+SAFE1 ) )
366            END IF
367   80    CONTINUE
368         BERR( J ) = S
369*
370*        Test stopping criterion. Continue iterating if
371*           1) The residual BERR(J) is larger than machine epsilon, and
372*           2) BERR(J) decreased by at least a factor of 2 during the
373*              last iteration, and
374*           3) At most ITMAX iterations tried.
375*
376         IF( BERR( J ).GT.EPS .AND. TWO*BERR( J ).LE.LSTRES .AND.
377     $       COUNT.LE.ITMAX ) THEN
378*
379*           Update solution and try again.
380*
381            CALL DGBTRS( TRANS, N, KL, KU, 1, AFB, LDAFB, IPIV,
382     $                   WORK( N+1 ), N, INFO )
383            CALL DAXPY( N, ONE, WORK( N+1 ), 1, X( 1, J ), 1 )
384            LSTRES = BERR( J )
385            COUNT = COUNT + 1
386            GO TO 20
387         END IF
388*
389*        Bound error from formula
390*
391*        norm(X - XTRUE) / norm(X) .le. FERR =
392*        norm( abs(inv(op(A)))*
393*           ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
394*
395*        where
396*          norm(Z) is the magnitude of the largest component of Z
397*          inv(op(A)) is the inverse of op(A)
398*          abs(Z) is the componentwise absolute value of the matrix or
399*             vector Z
400*          NZ is the maximum number of nonzeros in any row of A, plus 1
401*          EPS is machine epsilon
402*
403*        The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
404*        is incremented by SAFE1 if the i-th component of
405*        abs(op(A))*abs(X) + abs(B) is less than SAFE2.
406*
407*        Use DLACN2 to estimate the infinity-norm of the matrix
408*           inv(op(A)) * diag(W),
409*        where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
410*
411         DO 90 I = 1, N
412            IF( WORK( I ).GT.SAFE2 ) THEN
413               WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I )
414            ELSE
415               WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I ) + SAFE1
416            END IF
417   90    CONTINUE
418*
419         KASE = 0
420  100    CONTINUE
421         CALL DLACN2( N, WORK( 2*N+1 ), WORK( N+1 ), IWORK, FERR( J ),
422     $                KASE, ISAVE )
423         IF( KASE.NE.0 ) THEN
424            IF( KASE.EQ.1 ) THEN
425*
426*              Multiply by diag(W)*inv(op(A)**T).
427*
428               CALL DGBTRS( TRANST, N, KL, KU, 1, AFB, LDAFB, IPIV,
429     $                      WORK( N+1 ), N, INFO )
430               DO 110 I = 1, N
431                  WORK( N+I ) = WORK( N+I )*WORK( I )
432  110          CONTINUE
433            ELSE
434*
435*              Multiply by inv(op(A))*diag(W).
436*
437               DO 120 I = 1, N
438                  WORK( N+I ) = WORK( N+I )*WORK( I )
439  120          CONTINUE
440               CALL DGBTRS( TRANS, N, KL, KU, 1, AFB, LDAFB, IPIV,
441     $                      WORK( N+1 ), N, INFO )
442            END IF
443            GO TO 100
444         END IF
445*
446*        Normalize error.
447*
448         LSTRES = ZERO
449         DO 130 I = 1, N
450            LSTRES = MAX( LSTRES, ABS( X( I, J ) ) )
451  130    CONTINUE
452         IF( LSTRES.NE.ZERO )
453     $      FERR( J ) = FERR( J ) / LSTRES
454*
455  140 CONTINUE
456*
457      RETURN
458*
459*     End of DGBRFS
460*
461      END
462