```      SUBROUTINE DLAGTS( JOB, N, A, B, C, D, IN, Y, TOL, INFO )
*
*  -- LAPACK auxiliary routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
INTEGER            INFO, JOB, N
DOUBLE PRECISION   TOL
*     ..
*     .. Array Arguments ..
INTEGER            IN( * )
DOUBLE PRECISION   A( * ), B( * ), C( * ), D( * ), Y( * )
*     ..
*
*  Purpose
*  =======
*
*  DLAGTS may be used to solve one of the systems of equations
*
*     (T - lambda*I)*x = y   or   (T - lambda*I)'*x = y,
*
*  where T is an n by n tridiagonal matrix, for x, following the
*  factorization of (T - lambda*I) as
*
*     (T - lambda*I) = P*L*U ,
*
*  by routine DLAGTF. The choice of equation to be solved is
*  controlled by the argument JOB, and in each case there is an option
*  to perturb zero or very small diagonal elements of U, this option
*  being intended for use in applications such as inverse iteration.
*
*  Arguments
*  =========
*
*  JOB     (input) INTEGER
*          Specifies the job to be performed by DLAGTS as follows:
*          =  1: The equations  (T - lambda*I)x = y  are to be solved,
*                but diagonal elements of U are not to be perturbed.
*          = -1: The equations  (T - lambda*I)x = y  are to be solved
*                and, if overflow would otherwise occur, the diagonal
*                elements of U are to be perturbed. See argument TOL
*                below.
*          =  2: The equations  (T - lambda*I)'x = y  are to be solved,
*                but diagonal elements of U are not to be perturbed.
*          = -2: The equations  (T - lambda*I)'x = y  are to be solved
*                and, if overflow would otherwise occur, the diagonal
*                elements of U are to be perturbed. See argument TOL
*                below.
*
*  N       (input) INTEGER
*          The order of the matrix T.
*
*  A       (input) DOUBLE PRECISION array, dimension (N)
*          On entry, A must contain the diagonal elements of U as
*          returned from DLAGTF.
*
*  B       (input) DOUBLE PRECISION array, dimension (N-1)
*          On entry, B must contain the first super-diagonal elements of
*          U as returned from DLAGTF.
*
*  C       (input) DOUBLE PRECISION array, dimension (N-1)
*          On entry, C must contain the sub-diagonal elements of L as
*          returned from DLAGTF.
*
*  D       (input) DOUBLE PRECISION array, dimension (N-2)
*          On entry, D must contain the second super-diagonal elements
*          of U as returned from DLAGTF.
*
*  IN      (input) INTEGER array, dimension (N)
*          On entry, IN must contain details of the matrix P as returned
*          from DLAGTF.
*
*  Y       (input/output) DOUBLE PRECISION array, dimension (N)
*          On entry, the right hand side vector y.
*          On exit, Y is overwritten by the solution vector x.
*
*  TOL     (input/output) DOUBLE PRECISION
*          On entry, with  JOB .lt. 0, TOL should be the minimum
*          perturbation to be made to very small diagonal elements of U.
*          TOL should normally be chosen as about eps*norm(U), where eps
*          is the relative machine precision, but if TOL is supplied as
*          non-positive, then it is reset to eps*max( abs( u(i,j) ) ).
*          If  JOB .gt. 0  then TOL is not referenced.
*
*          On exit, TOL is changed as described above, only if TOL is
*          non-positive on entry. Otherwise TOL is unchanged.
*
*  INFO    (output) INTEGER
*          = 0   : successful exit
*          .lt. 0: if INFO = -i, the i-th argument had an illegal value
*          .gt. 0: overflow would occur when computing the INFO(th)
*                  element of the solution vector x. This can only occur
*                  when JOB is supplied as positive and either means
*                  that a diagonal element of U is very small, or that
*                  the elements of the right-hand side vector y are very
*                  large.
*
*  =====================================================================
*
*     .. Parameters ..
DOUBLE PRECISION   ONE, ZERO
PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
*     ..
*     .. Local Scalars ..
INTEGER            K
DOUBLE PRECISION   ABSAK, AK, BIGNUM, EPS, PERT, SFMIN, TEMP
*     ..
*     .. Intrinsic Functions ..
INTRINSIC          ABS, MAX, SIGN
*     ..
*     .. External Functions ..
DOUBLE PRECISION   DLAMCH
EXTERNAL           DLAMCH
*     ..
*     .. External Subroutines ..
EXTERNAL           XERBLA
*     ..
*     .. Executable Statements ..
*
INFO = 0
IF( ( ABS( JOB ).GT.2 ) .OR. ( JOB.EQ.0 ) ) THEN
INFO = -1
ELSE IF( N.LT.0 ) THEN
INFO = -2
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DLAGTS', -INFO )
RETURN
END IF
*
IF( N.EQ.0 )
\$   RETURN
*
EPS = DLAMCH( 'Epsilon' )
SFMIN = DLAMCH( 'Safe minimum' )
BIGNUM = ONE / SFMIN
*
IF( JOB.LT.0 ) THEN
IF( TOL.LE.ZERO ) THEN
TOL = ABS( A( 1 ) )
IF( N.GT.1 )
\$         TOL = MAX( TOL, ABS( A( 2 ) ), ABS( B( 1 ) ) )
DO 10 K = 3, N
TOL = MAX( TOL, ABS( A( K ) ), ABS( B( K-1 ) ),
\$               ABS( D( K-2 ) ) )
10       CONTINUE
TOL = TOL*EPS
IF( TOL.EQ.ZERO )
\$         TOL = EPS
END IF
END IF
*
IF( ABS( JOB ).EQ.1 ) THEN
DO 20 K = 2, N
IF( IN( K-1 ).EQ.0 ) THEN
Y( K ) = Y( K ) - C( K-1 )*Y( K-1 )
ELSE
TEMP = Y( K-1 )
Y( K-1 ) = Y( K )
Y( K ) = TEMP - C( K-1 )*Y( K )
END IF
20    CONTINUE
IF( JOB.EQ.1 ) THEN
DO 30 K = N, 1, -1
IF( K.LE.N-2 ) THEN
TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 )
ELSE IF( K.EQ.N-1 ) THEN
TEMP = Y( K ) - B( K )*Y( K+1 )
ELSE
TEMP = Y( K )
END IF
AK = A( K )
ABSAK = ABS( AK )
IF( ABSAK.LT.ONE ) THEN
IF( ABSAK.LT.SFMIN ) THEN
IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
\$                    THEN
INFO = K
RETURN
ELSE
TEMP = TEMP*BIGNUM
AK = AK*BIGNUM
END IF
ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
INFO = K
RETURN
END IF
END IF
Y( K ) = TEMP / AK
30       CONTINUE
ELSE
DO 50 K = N, 1, -1
IF( K.LE.N-2 ) THEN
TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 )
ELSE IF( K.EQ.N-1 ) THEN
TEMP = Y( K ) - B( K )*Y( K+1 )
ELSE
TEMP = Y( K )
END IF
AK = A( K )
PERT = SIGN( TOL, AK )
40          CONTINUE
ABSAK = ABS( AK )
IF( ABSAK.LT.ONE ) THEN
IF( ABSAK.LT.SFMIN ) THEN
IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
\$                    THEN
AK = AK + PERT
PERT = 2*PERT
GO TO 40
ELSE
TEMP = TEMP*BIGNUM
AK = AK*BIGNUM
END IF
ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
AK = AK + PERT
PERT = 2*PERT
GO TO 40
END IF
END IF
Y( K ) = TEMP / AK
50       CONTINUE
END IF
ELSE
*
*        Come to here if  JOB = 2 or -2
*
IF( JOB.EQ.2 ) THEN
DO 60 K = 1, N
IF( K.GE.3 ) THEN
TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 )
ELSE IF( K.EQ.2 ) THEN
TEMP = Y( K ) - B( K-1 )*Y( K-1 )
ELSE
TEMP = Y( K )
END IF
AK = A( K )
ABSAK = ABS( AK )
IF( ABSAK.LT.ONE ) THEN
IF( ABSAK.LT.SFMIN ) THEN
IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
\$                    THEN
INFO = K
RETURN
ELSE
TEMP = TEMP*BIGNUM
AK = AK*BIGNUM
END IF
ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
INFO = K
RETURN
END IF
END IF
Y( K ) = TEMP / AK
60       CONTINUE
ELSE
DO 80 K = 1, N
IF( K.GE.3 ) THEN
TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 )
ELSE IF( K.EQ.2 ) THEN
TEMP = Y( K ) - B( K-1 )*Y( K-1 )
ELSE
TEMP = Y( K )
END IF
AK = A( K )
PERT = SIGN( TOL, AK )
70          CONTINUE
ABSAK = ABS( AK )
IF( ABSAK.LT.ONE ) THEN
IF( ABSAK.LT.SFMIN ) THEN
IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
\$                    THEN
AK = AK + PERT
PERT = 2*PERT
GO TO 70
ELSE
TEMP = TEMP*BIGNUM
AK = AK*BIGNUM
END IF
ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
AK = AK + PERT
PERT = 2*PERT
GO TO 70
END IF
END IF
Y( K ) = TEMP / AK
80       CONTINUE
END IF
*
DO 90 K = N, 2, -1
IF( IN( K-1 ).EQ.0 ) THEN
Y( K-1 ) = Y( K-1 ) - C( K-1 )*Y( K )
ELSE
TEMP = Y( K-1 )
Y( K-1 ) = Y( K )
Y( K ) = TEMP - C( K-1 )*Y( K )
END IF
90    CONTINUE
END IF
*
*     End of DLAGTS
*
END

```