1*> \brief <b> DGTSVX computes the solution to system of linear equations A * X = B for GT matrices <b>
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DGTSVX + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgtsvx.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgtsvx.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgtsvx.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DGTSVX( FACT, TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF,
22*                          DU2, IPIV, B, LDB, X, LDX, RCOND, FERR, BERR,
23*                          WORK, IWORK, INFO )
24*
25*       .. Scalar Arguments ..
26*       CHARACTER          FACT, TRANS
27*       INTEGER            INFO, LDB, LDX, N, NRHS
28*       DOUBLE PRECISION   RCOND
29*       ..
30*       .. Array Arguments ..
31*       INTEGER            IPIV( * ), IWORK( * )
32*       DOUBLE PRECISION   B( LDB, * ), BERR( * ), D( * ), DF( * ),
33*      $                   DL( * ), DLF( * ), DU( * ), DU2( * ), DUF( * ),
34*      $                   FERR( * ), WORK( * ), X( LDX, * )
35*       ..
36*
37*
38*> \par Purpose:
39*  =============
40*>
41*> \verbatim
42*>
43*> DGTSVX uses the LU factorization to compute the solution to a real
44*> system of linear equations A * X = B or A**T * X = B,
45*> where A is a tridiagonal matrix of order N and X and B are N-by-NRHS
46*> matrices.
47*>
48*> Error bounds on the solution and a condition estimate are also
49*> provided.
50*> \endverbatim
51*
52*> \par Description:
53*  =================
54*>
55*> \verbatim
56*>
57*> The following steps are performed:
58*>
59*> 1. If FACT = 'N', the LU decomposition is used to factor the matrix A
60*>    as A = L * U, where L is a product of permutation and unit lower
61*>    bidiagonal matrices and U is upper triangular with nonzeros in
62*>    only the main diagonal and first two superdiagonals.
63*>
64*> 2. If some U(i,i)=0, so that U is exactly singular, then the routine
65*>    returns with INFO = i. Otherwise, the factored form of A is used
66*>    to estimate the condition number of the matrix A.  If the
67*>    reciprocal of the condition number is less than machine precision,
68*>    INFO = N+1 is returned as a warning, but the routine still goes on
69*>    to solve for X and compute error bounds as described below.
70*>
71*> 3. The system of equations is solved for X using the factored form
72*>    of A.
73*>
74*> 4. Iterative refinement is applied to improve the computed solution
75*>    matrix and calculate error bounds and backward error estimates
76*>    for it.
77*> \endverbatim
78*
79*  Arguments:
80*  ==========
81*
82*> \param[in] FACT
83*> \verbatim
84*>          FACT is CHARACTER*1
85*>          Specifies whether or not the factored form of A has been
86*>          supplied on entry.
87*>          = 'F':  DLF, DF, DUF, DU2, and IPIV contain the factored
88*>                  form of A; DL, D, DU, DLF, DF, DUF, DU2 and IPIV
89*>                  will not be modified.
90*>          = 'N':  The matrix will be copied to DLF, DF, and DUF
91*>                  and factored.
92*> \endverbatim
93*>
94*> \param[in] TRANS
95*> \verbatim
96*>          TRANS is CHARACTER*1
97*>          Specifies the form of the system of equations:
98*>          = 'N':  A * X = B     (No transpose)
99*>          = 'T':  A**T * X = B  (Transpose)
100*>          = 'C':  A**H * X = B  (Conjugate transpose = Transpose)
101*> \endverbatim
102*>
103*> \param[in] N
104*> \verbatim
105*>          N is INTEGER
106*>          The order of the matrix A.  N >= 0.
107*> \endverbatim
108*>
109*> \param[in] NRHS
110*> \verbatim
111*>          NRHS is INTEGER
112*>          The number of right hand sides, i.e., the number of columns
113*>          of the matrix B.  NRHS >= 0.
114*> \endverbatim
115*>
116*> \param[in] DL
117*> \verbatim
118*>          DL is DOUBLE PRECISION array, dimension (N-1)
119*>          The (n-1) subdiagonal elements of A.
120*> \endverbatim
121*>
122*> \param[in] D
123*> \verbatim
124*>          D is DOUBLE PRECISION array, dimension (N)
125*>          The n diagonal elements of A.
126*> \endverbatim
127*>
128*> \param[in] DU
129*> \verbatim
130*>          DU is DOUBLE PRECISION array, dimension (N-1)
131*>          The (n-1) superdiagonal elements of A.
132*> \endverbatim
133*>
134*> \param[in,out] DLF
135*> \verbatim
136*>          DLF is DOUBLE PRECISION array, dimension (N-1)
137*>          If FACT = 'F', then DLF is an input argument and on entry
138*>          contains the (n-1) multipliers that define the matrix L from
139*>          the LU factorization of A as computed by DGTTRF.
140*>
141*>          If FACT = 'N', then DLF is an output argument and on exit
142*>          contains the (n-1) multipliers that define the matrix L from
143*>          the LU factorization of A.
144*> \endverbatim
145*>
146*> \param[in,out] DF
147*> \verbatim
148*>          DF is DOUBLE PRECISION array, dimension (N)
149*>          If FACT = 'F', then DF is an input argument and on entry
150*>          contains the n diagonal elements of the upper triangular
151*>          matrix U from the LU factorization of A.
152*>
153*>          If FACT = 'N', then DF is an output argument and on exit
154*>          contains the n diagonal elements of the upper triangular
155*>          matrix U from the LU factorization of A.
156*> \endverbatim
157*>
158*> \param[in,out] DUF
159*> \verbatim
160*>          DUF is DOUBLE PRECISION array, dimension (N-1)
161*>          If FACT = 'F', then DUF is an input argument and on entry
162*>          contains the (n-1) elements of the first superdiagonal of U.
163*>
164*>          If FACT = 'N', then DUF is an output argument and on exit
165*>          contains the (n-1) elements of the first superdiagonal of U.
166*> \endverbatim
167*>
168*> \param[in,out] DU2
169*> \verbatim
170*>          DU2 is DOUBLE PRECISION array, dimension (N-2)
171*>          If FACT = 'F', then DU2 is an input argument and on entry
172*>          contains the (n-2) elements of the second superdiagonal of
173*>          U.
174*>
175*>          If FACT = 'N', then DU2 is an output argument and on exit
176*>          contains the (n-2) elements of the second superdiagonal of
177*>          U.
178*> \endverbatim
179*>
180*> \param[in,out] IPIV
181*> \verbatim
182*>          IPIV is INTEGER array, dimension (N)
183*>          If FACT = 'F', then IPIV is an input argument and on entry
184*>          contains the pivot indices from the LU factorization of A as
185*>          computed by DGTTRF.
186*>
187*>          If FACT = 'N', then IPIV is an output argument and on exit
188*>          contains the pivot indices from the LU factorization of A;
189*>          row i of the matrix was interchanged with row IPIV(i).
190*>          IPIV(i) will always be either i or i+1; IPIV(i) = i indicates
191*>          a row interchange was not required.
192*> \endverbatim
193*>
194*> \param[in] B
195*> \verbatim
196*>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
197*>          The N-by-NRHS right hand side matrix B.
198*> \endverbatim
199*>
200*> \param[in] LDB
201*> \verbatim
202*>          LDB is INTEGER
203*>          The leading dimension of the array B.  LDB >= max(1,N).
204*> \endverbatim
205*>
206*> \param[out] X
207*> \verbatim
208*>          X is DOUBLE PRECISION array, dimension (LDX,NRHS)
209*>          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
210*> \endverbatim
211*>
212*> \param[in] LDX
213*> \verbatim
214*>          LDX is INTEGER
215*>          The leading dimension of the array X.  LDX >= max(1,N).
216*> \endverbatim
217*>
218*> \param[out] RCOND
219*> \verbatim
220*>          RCOND is DOUBLE PRECISION
221*>          The estimate of the reciprocal condition number of the matrix
222*>          A.  If RCOND is less than the machine precision (in
223*>          particular, if RCOND = 0), the matrix is singular to working
224*>          precision.  This condition is indicated by a return code of
225*>          INFO > 0.
226*> \endverbatim
227*>
228*> \param[out] FERR
229*> \verbatim
230*>          FERR is DOUBLE PRECISION array, dimension (NRHS)
231*>          The estimated forward error bound for each solution vector
232*>          X(j) (the j-th column of the solution matrix X).
233*>          If XTRUE is the true solution corresponding to X(j), FERR(j)
234*>          is an estimated upper bound for the magnitude of the largest
235*>          element in (X(j) - XTRUE) divided by the magnitude of the
236*>          largest element in X(j).  The estimate is as reliable as
237*>          the estimate for RCOND, and is almost always a slight
238*>          overestimate of the true error.
239*> \endverbatim
240*>
241*> \param[out] BERR
242*> \verbatim
243*>          BERR is DOUBLE PRECISION array, dimension (NRHS)
244*>          The componentwise relative backward error of each solution
245*>          vector X(j) (i.e., the smallest relative change in
246*>          any element of A or B that makes X(j) an exact solution).
247*> \endverbatim
248*>
249*> \param[out] WORK
250*> \verbatim
251*>          WORK is DOUBLE PRECISION array, dimension (3*N)
252*> \endverbatim
253*>
254*> \param[out] IWORK
255*> \verbatim
256*>          IWORK is INTEGER array, dimension (N)
257*> \endverbatim
258*>
259*> \param[out] INFO
260*> \verbatim
261*>          INFO is INTEGER
262*>          = 0:  successful exit
263*>          < 0:  if INFO = -i, the i-th argument had an illegal value
264*>          > 0:  if INFO = i, and i is
265*>                <= N:  U(i,i) is exactly zero.  The factorization
266*>                       has not been completed unless i = N, but the
267*>                       factor U is exactly singular, so the solution
268*>                       and error bounds could not be computed.
269*>                       RCOND = 0 is returned.
270*>                = N+1: U is nonsingular, but RCOND is less than machine
271*>                       precision, meaning that the matrix is singular
272*>                       to working precision.  Nevertheless, the
273*>                       solution and error bounds are computed because
274*>                       there are a number of situations where the
275*>                       computed solution can be more accurate than the
276*>                       value of RCOND would suggest.
277*> \endverbatim
278*
279*  Authors:
280*  ========
281*
282*> \author Univ. of Tennessee
283*> \author Univ. of California Berkeley
284*> \author Univ. of Colorado Denver
285*> \author NAG Ltd.
286*
287*> \date September 2012
288*
289*> \ingroup doubleGTsolve
290*
291*  =====================================================================
292      SUBROUTINE DGTSVX( FACT, TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF,
293     $                   DU2, IPIV, B, LDB, X, LDX, RCOND, FERR, BERR,
294     $                   WORK, IWORK, INFO )
295*
296*  -- LAPACK driver routine (version 3.4.2) --
297*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
298*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
299*     September 2012
300*
301*     .. Scalar Arguments ..
302      CHARACTER          FACT, TRANS
303      INTEGER            INFO, LDB, LDX, N, NRHS
304      DOUBLE PRECISION   RCOND
305*     ..
306*     .. Array Arguments ..
307      INTEGER            IPIV( * ), IWORK( * )
308      DOUBLE PRECISION   B( LDB, * ), BERR( * ), D( * ), DF( * ),
309     $                   DL( * ), DLF( * ), DU( * ), DU2( * ), DUF( * ),
310     $                   FERR( * ), WORK( * ), X( LDX, * )
311*     ..
312*
313*  =====================================================================
314*
315*     .. Parameters ..
316      DOUBLE PRECISION   ZERO
317      PARAMETER          ( ZERO = 0.0D+0 )
318*     ..
319*     .. Local Scalars ..
320      LOGICAL            NOFACT, NOTRAN
321      CHARACTER          NORM
322      DOUBLE PRECISION   ANORM
323*     ..
324*     .. External Functions ..
325      LOGICAL            LSAME
326      DOUBLE PRECISION   DLAMCH, DLANGT
327      EXTERNAL           LSAME, DLAMCH, DLANGT
328*     ..
329*     .. External Subroutines ..
330      EXTERNAL           DCOPY, DGTCON, DGTRFS, DGTTRF, DGTTRS, DLACPY,
331     $                   XERBLA
332*     ..
333*     .. Intrinsic Functions ..
334      INTRINSIC          MAX
335*     ..
336*     .. Executable Statements ..
337*
338      INFO = 0
339      NOFACT = LSAME( FACT, 'N' )
340      NOTRAN = LSAME( TRANS, 'N' )
341      IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
342         INFO = -1
343      ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
344     $         LSAME( TRANS, 'C' ) ) THEN
345         INFO = -2
346      ELSE IF( N.LT.0 ) THEN
347         INFO = -3
348      ELSE IF( NRHS.LT.0 ) THEN
349         INFO = -4
350      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
351         INFO = -14
352      ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
353         INFO = -16
354      END IF
355      IF( INFO.NE.0 ) THEN
356         CALL XERBLA( 'DGTSVX', -INFO )
357         RETURN
358      END IF
359*
360      IF( NOFACT ) THEN
361*
362*        Compute the LU factorization of A.
363*
364         CALL DCOPY( N, D, 1, DF, 1 )
365         IF( N.GT.1 ) THEN
366            CALL DCOPY( N-1, DL, 1, DLF, 1 )
367            CALL DCOPY( N-1, DU, 1, DUF, 1 )
368         END IF
369         CALL DGTTRF( N, DLF, DF, DUF, DU2, IPIV, INFO )
370*
371*        Return if INFO is non-zero.
372*
373         IF( INFO.GT.0 )THEN
374            RCOND = ZERO
375            RETURN
376         END IF
377      END IF
378*
379*     Compute the norm of the matrix A.
380*
381      IF( NOTRAN ) THEN
382         NORM = '1'
383      ELSE
384         NORM = 'I'
385      END IF
386      ANORM = DLANGT( NORM, N, DL, D, DU )
387*
388*     Compute the reciprocal of the condition number of A.
389*
390      CALL DGTCON( NORM, N, DLF, DF, DUF, DU2, IPIV, ANORM, RCOND, WORK,
391     $             IWORK, INFO )
392*
393*     Compute the solution vectors X.
394*
395      CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
396      CALL DGTTRS( TRANS, N, NRHS, DLF, DF, DUF, DU2, IPIV, X, LDX,
397     $             INFO )
398*
399*     Use iterative refinement to improve the computed solutions and
400*     compute error bounds and backward error estimates for them.
401*
402      CALL DGTRFS( TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF, DU2, IPIV,
403     $             B, LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO )
404*
405*     Set INFO = N+1 if the matrix is singular to working precision.
406*
407      IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
408     $   INFO = N + 1
409*
410      RETURN
411*
412*     End of DGTSVX
413*
414      END
415