1*> \brief \b DLAGTF computes an LU factorization of a matrix T-λI, where T is a general tridiagonal matrix, and λ a scalar, using partial pivoting with row interchanges.
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DLAGTF + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlagtf.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlagtf.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlagtf.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DLAGTF( N, A, LAMBDA, B, C, TOL, D, IN, INFO )
22*
23*       .. Scalar Arguments ..
24*       INTEGER            INFO, N
25*       DOUBLE PRECISION   LAMBDA, TOL
26*       ..
27*       .. Array Arguments ..
28*       INTEGER            IN( * )
29*       DOUBLE PRECISION   A( * ), B( * ), C( * ), D( * )
30*       ..
31*
32*
33*> \par Purpose:
34*  =============
35*>
36*> \verbatim
37*>
38*> DLAGTF factorizes the matrix (T - lambda*I), where T is an n by n
39*> tridiagonal matrix and lambda is a scalar, as
40*>
41*>    T - lambda*I = PLU,
42*>
43*> where P is a permutation matrix, L is a unit lower tridiagonal matrix
44*> with at most one non-zero sub-diagonal elements per column and U is
45*> an upper triangular matrix with at most two non-zero super-diagonal
46*> elements per column.
47*>
48*> The factorization is obtained by Gaussian elimination with partial
49*> pivoting and implicit row scaling.
50*>
51*> The parameter LAMBDA is included in the routine so that DLAGTF may
52*> be used, in conjunction with DLAGTS, to obtain eigenvectors of T by
53*> inverse iteration.
54*> \endverbatim
55*
56*  Arguments:
57*  ==========
58*
59*> \param[in] N
60*> \verbatim
61*>          N is INTEGER
62*>          The order of the matrix T.
63*> \endverbatim
64*>
65*> \param[in,out] A
66*> \verbatim
67*>          A is DOUBLE PRECISION array, dimension (N)
68*>          On entry, A must contain the diagonal elements of T.
69*>
70*>          On exit, A is overwritten by the n diagonal elements of the
71*>          upper triangular matrix U of the factorization of T.
72*> \endverbatim
73*>
74*> \param[in] LAMBDA
75*> \verbatim
76*>          LAMBDA is DOUBLE PRECISION
77*>          On entry, the scalar lambda.
78*> \endverbatim
79*>
80*> \param[in,out] B
81*> \verbatim
82*>          B is DOUBLE PRECISION array, dimension (N-1)
83*>          On entry, B must contain the (n-1) super-diagonal elements of
84*>          T.
85*>
86*>          On exit, B is overwritten by the (n-1) super-diagonal
87*>          elements of the matrix U of the factorization of T.
88*> \endverbatim
89*>
90*> \param[in,out] C
91*> \verbatim
92*>          C is DOUBLE PRECISION array, dimension (N-1)
93*>          On entry, C must contain the (n-1) sub-diagonal elements of
94*>          T.
95*>
96*>          On exit, C is overwritten by the (n-1) sub-diagonal elements
97*>          of the matrix L of the factorization of T.
98*> \endverbatim
99*>
100*> \param[in] TOL
101*> \verbatim
102*>          TOL is DOUBLE PRECISION
103*>          On entry, a relative tolerance used to indicate whether or
104*>          not the matrix (T - lambda*I) is nearly singular. TOL should
105*>          normally be chose as approximately the largest relative error
106*>          in the elements of T. For example, if the elements of T are
107*>          correct to about 4 significant figures, then TOL should be
108*>          set to about 5*10**(-4). If TOL is supplied as less than eps,
109*>          where eps is the relative machine precision, then the value
110*>          eps is used in place of TOL.
111*> \endverbatim
112*>
113*> \param[out] D
114*> \verbatim
115*>          D is DOUBLE PRECISION array, dimension (N-2)
116*>          On exit, D is overwritten by the (n-2) second super-diagonal
117*>          elements of the matrix U of the factorization of T.
118*> \endverbatim
119*>
120*> \param[out] IN
121*> \verbatim
122*>          IN is INTEGER array, dimension (N)
123*>          On exit, IN contains details of the permutation matrix P. If
124*>          an interchange occurred at the kth step of the elimination,
125*>          then IN(k) = 1, otherwise IN(k) = 0. The element IN(n)
126*>          returns the smallest positive integer j such that
127*>
128*>             abs( u(j,j) ) <= norm( (T - lambda*I)(j) )*TOL,
129*>
130*>          where norm( A(j) ) denotes the sum of the absolute values of
131*>          the jth row of the matrix A. If no such j exists then IN(n)
132*>          is returned as zero. If IN(n) is returned as positive, then a
133*>          diagonal element of U is small, indicating that
134*>          (T - lambda*I) is singular or nearly singular,
135*> \endverbatim
136*>
137*> \param[out] INFO
138*> \verbatim
139*>          INFO is INTEGER
140*>          = 0:  successful exit
141*>          < 0:  if INFO = -k, the kth argument had an illegal value
142*> \endverbatim
143*
144*  Authors:
145*  ========
146*
147*> \author Univ. of Tennessee
148*> \author Univ. of California Berkeley
149*> \author Univ. of Colorado Denver
150*> \author NAG Ltd.
151*
152*> \date December 2016
153*
154*> \ingroup auxOTHERcomputational
155*
156*  =====================================================================
157      SUBROUTINE DLAGTF( N, A, LAMBDA, B, C, TOL, D, IN, INFO )
158*
159*  -- LAPACK computational routine (version 3.7.0) --
160*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
161*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
162*     December 2016
163*
164*     .. Scalar Arguments ..
165      INTEGER            INFO, N
166      DOUBLE PRECISION   LAMBDA, TOL
167*     ..
168*     .. Array Arguments ..
169      INTEGER            IN( * )
170      DOUBLE PRECISION   A( * ), B( * ), C( * ), D( * )
171*     ..
172*
173* =====================================================================
174*
175*     .. Parameters ..
176      DOUBLE PRECISION   ZERO
177      PARAMETER          ( ZERO = 0.0D+0 )
178*     ..
179*     .. Local Scalars ..
180      INTEGER            K
181      DOUBLE PRECISION   EPS, MULT, PIV1, PIV2, SCALE1, SCALE2, TEMP, TL
182*     ..
183*     .. Intrinsic Functions ..
184      INTRINSIC          ABS, MAX
185*     ..
186*     .. External Functions ..
187      DOUBLE PRECISION   DLAMCH
188      EXTERNAL           DLAMCH
189*     ..
190*     .. External Subroutines ..
191      EXTERNAL           XERBLA
192*     ..
193*     .. Executable Statements ..
194*
195      INFO = 0
196      IF( N.LT.0 ) THEN
197         INFO = -1
198         CALL XERBLA( 'DLAGTF', -INFO )
199         RETURN
200      END IF
201*
202      IF( N.EQ.0 )
203     $   RETURN
204*
205      A( 1 ) = A( 1 ) - LAMBDA
206      IN( N ) = 0
207      IF( N.EQ.1 ) THEN
208         IF( A( 1 ).EQ.ZERO )
209     $      IN( 1 ) = 1
210         RETURN
211      END IF
212*
213      EPS = DLAMCH( 'Epsilon' )
214*
215      TL = MAX( TOL, EPS )
216      SCALE1 = ABS( A( 1 ) ) + ABS( B( 1 ) )
217      DO 10 K = 1, N - 1
218         A( K+1 ) = A( K+1 ) - LAMBDA
219         SCALE2 = ABS( C( K ) ) + ABS( A( K+1 ) )
220         IF( K.LT.( N-1 ) )
221     $      SCALE2 = SCALE2 + ABS( B( K+1 ) )
222         IF( A( K ).EQ.ZERO ) THEN
223            PIV1 = ZERO
224         ELSE
225            PIV1 = ABS( A( K ) ) / SCALE1
226         END IF
227         IF( C( K ).EQ.ZERO ) THEN
228            IN( K ) = 0
229            PIV2 = ZERO
230            SCALE1 = SCALE2
231            IF( K.LT.( N-1 ) )
232     $         D( K ) = ZERO
233         ELSE
234            PIV2 = ABS( C( K ) ) / SCALE2
235            IF( PIV2.LE.PIV1 ) THEN
236               IN( K ) = 0
237               SCALE1 = SCALE2
238               C( K ) = C( K ) / A( K )
239               A( K+1 ) = A( K+1 ) - C( K )*B( K )
240               IF( K.LT.( N-1 ) )
241     $            D( K ) = ZERO
242            ELSE
243               IN( K ) = 1
244               MULT = A( K ) / C( K )
245               A( K ) = C( K )
246               TEMP = A( K+1 )
247               A( K+1 ) = B( K ) - MULT*TEMP
248               IF( K.LT.( N-1 ) ) THEN
249                  D( K ) = B( K+1 )
250                  B( K+1 ) = -MULT*D( K )
251               END IF
252               B( K ) = TEMP
253               C( K ) = MULT
254            END IF
255         END IF
256         IF( ( MAX( PIV1, PIV2 ).LE.TL ) .AND. ( IN( N ).EQ.0 ) )
257     $      IN( N ) = K
258   10 CONTINUE
259      IF( ( ABS( A( N ) ).LE.SCALE1*TL ) .AND. ( IN( N ).EQ.0 ) )
260     $   IN( N ) = N
261*
262      RETURN
263*
264*     End of DLAGTF
265*
266      END
267