dlagts.f

Go to the documentation of this file.

*> \brief \b DLAGTS solves the system of equations (T-λI)x = y or (T-λI)Tx = y,where T is a general tridiagonal matrix and λ a scalar, using the LU factorization computed by slagtf.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*> \htmlonly
*> Download DLAGTS + dependencies 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlagts.f"> 
*> [TGZ]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlagts.f"> 
*> [ZIP]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlagts.f"> 
*> [TXT]</a>
*> \endhtmlonly 
*
*  Definition:
*  ===========
*
*       SUBROUTINE DLAGTS( JOB, N, A, B, C, D, IN, Y, TOL, INFO )
* 
*       .. Scalar Arguments ..
*       INTEGER            INFO, JOB, N
*       DOUBLE PRECISION   TOL
*       ..
*       .. Array Arguments ..
*       INTEGER            IN( * )
*       DOUBLE PRECISION   A( * ), B( * ), C( * ), D( * ), Y( * )
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DLAGTS may be used to solve one of the systems of equations
*>
*>    (T - lambda*I)*x = y   or   (T - lambda*I)**T*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.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] JOB
*> \verbatim
*>          JOB is 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)**Tx = y  are to be solved,
*>                but diagonal elements of U are not to be perturbed.
*>          = -2: The equations  (T - lambda*I)**Tx = y  are to be solved
*>                and, if overflow would otherwise occur, the diagonal
*>                elements of U are to be perturbed. See argument TOL
*>                below.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrix T.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*>          A is DOUBLE PRECISION array, dimension (N)
*>          On entry, A must contain the diagonal elements of U as
*>          returned from DLAGTF.
*> \endverbatim
*>
*> \param[in] B
*> \verbatim
*>          B is DOUBLE PRECISION array, dimension (N-1)
*>          On entry, B must contain the first super-diagonal elements of
*>          U as returned from DLAGTF.
*> \endverbatim
*>
*> \param[in] C
*> \verbatim
*>          C is DOUBLE PRECISION array, dimension (N-1)
*>          On entry, C must contain the sub-diagonal elements of L as
*>          returned from DLAGTF.
*> \endverbatim
*>
*> \param[in] D
*> \verbatim
*>          D is DOUBLE PRECISION array, dimension (N-2)
*>          On entry, D must contain the second super-diagonal elements
*>          of U as returned from DLAGTF.
*> \endverbatim
*>
*> \param[in] IN
*> \verbatim
*>          IN is INTEGER array, dimension (N)
*>          On entry, IN must contain details of the matrix P as returned
*>          from DLAGTF.
*> \endverbatim
*>
*> \param[in,out] Y
*> \verbatim
*>          Y is DOUBLE PRECISION array, dimension (N)
*>          On entry, the right hand side vector y.
*>          On exit, Y is overwritten by the solution vector x.
*> \endverbatim
*>
*> \param[in,out] TOL
*> \verbatim
*>          TOL is 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.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is 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.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee 
*> \author Univ. of California Berkeley 
*> \author Univ. of Colorado Denver 
*> \author NAG Ltd. 
*
*> \date September 2012
*
*> \ingroup auxOTHERauxiliary
*
*  =====================================================================
      SUBROUTINE dlagts( JOB, N, A, B, C, D, IN, Y, TOL, INFO )
*
*  -- LAPACK auxiliary routine (version 3.4.2) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     September 2012
*
*     .. Scalar Arguments ..
      INTEGER            INFO, JOB, N
      DOUBLE PRECISION   TOL
*     ..
*     .. Array Arguments ..
      INTEGER            IN( * )
      DOUBLE PRECISION   A( * ), B( * ), C( * ), D( * ), Y( * )
*     ..
*
*  =====================================================================
*
*     .. 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