1*> \brief <b> SGTSVX 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 SGTSVX + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgtsvx.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgtsvx.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgtsvx.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE SGTSVX( 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*       REAL               RCOND
29*       ..
30*       .. Array Arguments ..
31*       INTEGER            IPIV( * ), IWORK( * )
32*       REAL               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*> SGTSVX 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 REAL array, dimension (N-1)
119*>          The (n-1) subdiagonal elements of A.
120*> \endverbatim
121*>
122*> \param[in] D
123*> \verbatim
124*>          D is REAL array, dimension (N)
125*>          The n diagonal elements of A.
126*> \endverbatim
127*>
128*> \param[in] DU
129*> \verbatim
130*>          DU is REAL 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 REAL 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 SGTTRF.
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 REAL 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 REAL 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 REAL 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 SGTTRF.
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 REAL 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 REAL 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 REAL
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 REAL 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 REAL 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 REAL 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*> \ingroup realGTsolve
288*
289*  =====================================================================
290      SUBROUTINE SGTSVX( FACT, TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF,
291     $                   DU2, IPIV, B, LDB, X, LDX, RCOND, FERR, BERR,
292     $                   WORK, IWORK, INFO )
293*
294*  -- LAPACK driver routine --
295*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
296*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
297*
298*     .. Scalar Arguments ..
299      CHARACTER          FACT, TRANS
300      INTEGER            INFO, LDB, LDX, N, NRHS
301      REAL               RCOND
302*     ..
303*     .. Array Arguments ..
304      INTEGER            IPIV( * ), IWORK( * )
305      REAL               B( LDB, * ), BERR( * ), D( * ), DF( * ),
306     $                   DL( * ), DLF( * ), DU( * ), DU2( * ), DUF( * ),
307     $                   FERR( * ), WORK( * ), X( LDX, * )
308*     ..
309*
310*  =====================================================================
311*
312*     .. Parameters ..
313      REAL               ZERO
314      PARAMETER          ( ZERO = 0.0E+0 )
315*     ..
316*     .. Local Scalars ..
317      LOGICAL            NOFACT, NOTRAN
318      CHARACTER          NORM
319      REAL               ANORM
320*     ..
321*     .. External Functions ..
322      LOGICAL            LSAME
323      REAL               SLAMCH, SLANGT
324      EXTERNAL           LSAME, SLAMCH, SLANGT
325*     ..
326*     .. External Subroutines ..
327      EXTERNAL           SCOPY, SGTCON, SGTRFS, SGTTRF, SGTTRS, SLACPY,
328     $                   XERBLA
329*     ..
330*     .. Intrinsic Functions ..
331      INTRINSIC          MAX
332*     ..
333*     .. Executable Statements ..
334*
335      INFO = 0
336      NOFACT = LSAME( FACT, 'N' )
337      NOTRAN = LSAME( TRANS, 'N' )
338      IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
339         INFO = -1
340      ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
341     $         LSAME( TRANS, 'C' ) ) THEN
342         INFO = -2
343      ELSE IF( N.LT.0 ) THEN
344         INFO = -3
345      ELSE IF( NRHS.LT.0 ) THEN
346         INFO = -4
347      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
348         INFO = -14
349      ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
350         INFO = -16
351      END IF
352      IF( INFO.NE.0 ) THEN
353         CALL XERBLA( 'SGTSVX', -INFO )
354         RETURN
355      END IF
356*
357      IF( NOFACT ) THEN
358*
359*        Compute the LU factorization of A.
360*
361         CALL SCOPY( N, D, 1, DF, 1 )
362         IF( N.GT.1 ) THEN
363            CALL SCOPY( N-1, DL, 1, DLF, 1 )
364            CALL SCOPY( N-1, DU, 1, DUF, 1 )
365         END IF
366         CALL SGTTRF( N, DLF, DF, DUF, DU2, IPIV, INFO )
367*
368*        Return if INFO is non-zero.
369*
370         IF( INFO.GT.0 )THEN
371            RCOND = ZERO
372            RETURN
373         END IF
374      END IF
375*
376*     Compute the norm of the matrix A.
377*
378      IF( NOTRAN ) THEN
379         NORM = '1'
380      ELSE
381         NORM = 'I'
382      END IF
383      ANORM = SLANGT( NORM, N, DL, D, DU )
384*
385*     Compute the reciprocal of the condition number of A.
386*
387      CALL SGTCON( NORM, N, DLF, DF, DUF, DU2, IPIV, ANORM, RCOND, WORK,
388     $             IWORK, INFO )
389*
390*     Compute the solution vectors X.
391*
392      CALL SLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
393      CALL SGTTRS( TRANS, N, NRHS, DLF, DF, DUF, DU2, IPIV, X, LDX,
394     $             INFO )
395*
396*     Use iterative refinement to improve the computed solutions and
397*     compute error bounds and backward error estimates for them.
398*
399      CALL SGTRFS( TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF, DU2, IPIV,
400     $             B, LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO )
401*
402*     Set INFO = N+1 if the matrix is singular to working precision.
403*
404      IF( RCOND.LT.SLAMCH( 'Epsilon' ) )
405     $   INFO = N + 1
406*
407      RETURN
408*
409*     End of SGTSVX
410*
411      END
412