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