dlagtf.f

Go to the documentation of this file.

*> \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.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DLAGTF + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlagtf.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlagtf.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlagtf.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE DLAGTF( N, A, LAMBDA, B, C, TOL, D, IN, INFO )
*
*       .. Scalar Arguments ..
*       INTEGER            INFO, N
*       DOUBLE PRECISION   LAMBDA, TOL
*       ..
*       .. Array Arguments ..
*       INTEGER            IN( * )
*       DOUBLE PRECISION   A( * ), B( * ), C( * ), D( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DLAGTF factorizes the matrix (T - lambda*I), where T is an n by n
*> tridiagonal matrix and lambda is a scalar, as
*>
*>    T - lambda*I = PLU,
*>
*> where P is a permutation matrix, L is a unit lower tridiagonal matrix
*> with at most one non-zero sub-diagonal elements per column and U is
*> an upper triangular matrix with at most two non-zero super-diagonal
*> elements per column.
*>
*> The factorization is obtained by Gaussian elimination with partial
*> pivoting and implicit row scaling.
*>
*> The parameter LAMBDA is included in the routine so that DLAGTF may
*> be used, in conjunction with DLAGTS, to obtain eigenvectors of T by
*> inverse iteration.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrix T.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*>          A is DOUBLE PRECISION array, dimension (N)
*>          On entry, A must contain the diagonal elements of T.
*>
*>          On exit, A is overwritten by the n diagonal elements of the
*>          upper triangular matrix U of the factorization of T.
*> \endverbatim
*>
*> \param[in] LAMBDA
*> \verbatim
*>          LAMBDA is DOUBLE PRECISION
*>          On entry, the scalar lambda.
*> \endverbatim
*>
*> \param[in,out] B
*> \verbatim
*>          B is DOUBLE PRECISION array, dimension (N-1)
*>          On entry, B must contain the (n-1) super-diagonal elements of
*>          T.
*>
*>          On exit, B is overwritten by the (n-1) super-diagonal
*>          elements of the matrix U of the factorization of T.
*> \endverbatim
*>
*> \param[in,out] C
*> \verbatim
*>          C is DOUBLE PRECISION array, dimension (N-1)
*>          On entry, C must contain the (n-1) sub-diagonal elements of
*>          T.
*>
*>          On exit, C is overwritten by the (n-1) sub-diagonal elements
*>          of the matrix L of the factorization of T.
*> \endverbatim
*>
*> \param[in] TOL
*> \verbatim
*>          TOL is DOUBLE PRECISION
*>          On entry, a relative tolerance used to indicate whether or
*>          not the matrix (T - lambda*I) is nearly singular. TOL should
*>          normally be chose as approximately the largest relative error
*>          in the elements of T. For example, if the elements of T are
*>          correct to about 4 significant figures, then TOL should be
*>          set to about 5*10**(-4). If TOL is supplied as less than eps,
*>          where eps is the relative machine precision, then the value
*>          eps is used in place of TOL.
*> \endverbatim
*>
*> \param[out] D
*> \verbatim
*>          D is DOUBLE PRECISION array, dimension (N-2)
*>          On exit, D is overwritten by the (n-2) second super-diagonal
*>          elements of the matrix U of the factorization of T.
*> \endverbatim
*>
*> \param[out] IN
*> \verbatim
*>          IN is INTEGER array, dimension (N)
*>          On exit, IN contains details of the permutation matrix P. If
*>          an interchange occurred at the kth step of the elimination,
*>          then IN(k) = 1, otherwise IN(k) = 0. The element IN(n)
*>          returns the smallest positive integer j such that
*>
*>             abs( u(j,j) ) <= norm( (T - lambda*I)(j) )*TOL,
*>
*>          where norm( A(j) ) denotes the sum of the absolute values of
*>          the jth row of the matrix A. If no such j exists then IN(n)
*>          is returned as zero. If IN(n) is returned as positive, then a
*>          diagonal element of U is small, indicating that
*>          (T - lambda*I) is singular or nearly singular,
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0:  successful exit
*>          < 0:  if INFO = -k, the kth argument had an illegal value
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup lagtf
*
*  =====================================================================
      SUBROUTINE dlagtf( N, A, LAMBDA, B, C, TOL, D, IN, INFO )
*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER            INFO, N
      DOUBLE PRECISION   LAMBDA, TOL
*     ..
*     .. Array Arguments ..
      INTEGER            IN( * )
      DOUBLE PRECISION   A( * ), B( * ), C( * ), D( * )
*     ..
*
* =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO
      parameter( zero = 0.0d+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            K
      DOUBLE PRECISION   EPS, MULT, PIV1, PIV2, SCALE1, SCALE2, TEMP, TL
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH
      EXTERNAL           dlamch
*     ..
*     .. External Subroutines ..
      EXTERNAL           xerbla
*     ..
*     .. Executable Statements ..
*
      info = 0
      IF( n.LT.0 ) THEN
         info = -1
         CALL xerbla( 'DLAGTF', -info )
         RETURN
      END IF
*
      IF( n.EQ.0 )
     $   RETURN
*
      a( 1 ) = a( 1 ) - lambda
      in( n ) = 0
      IF( n.EQ.1 ) THEN
         IF( a( 1 ).EQ.zero )
     $      in( 1 ) = 1
         RETURN
      END IF
*
      eps = dlamch( 'Epsilon' )
*
      tl = max( tol, eps )
      scale1 = abs( a( 1 ) ) + abs( b( 1 ) )
      DO 10 k = 1, n - 1
         a( k+1 ) = a( k+1 ) - lambda
         scale2 = abs( c( k ) ) + abs( a( k+1 ) )
         IF( k.LT.( n-1 ) )
     $      scale2 = scale2 + abs( b( k+1 ) )
         IF( a( k ).EQ.zero ) THEN
            piv1 = zero
         ELSE
            piv1 = abs( a( k ) ) / scale1
         END IF
         IF( c( k ).EQ.zero ) THEN
            in( k ) = 0
            piv2 = zero
            scale1 = scale2
            IF( k.LT.( n-1 ) )
     $         d( k ) = zero
         ELSE
            piv2 = abs( c( k ) ) / scale2
            IF( piv2.LE.piv1 ) THEN
               in( k ) = 0
               scale1 = scale2
               c( k ) = c( k ) / a( k )
               a( k+1 ) = a( k+1 ) - c( k )*b( k )
               IF( k.LT.( n-1 ) )
     $            d( k ) = zero
            ELSE
               in( k ) = 1
               mult = a( k ) / c( k )
               a( k ) = c( k )
               temp = a( k+1 )
               a( k+1 ) = b( k ) - mult*temp
               IF( k.LT.( n-1 ) ) THEN
                  d( k ) = b( k+1 )
                  b( k+1 ) = -mult*d( k )
               END IF
               b( k ) = temp
               c( k ) = mult
            END IF
         END IF
         IF( ( max( piv1, piv2 ).LE.tl ) .AND. ( in( n ).EQ.0 ) )
     $      in( n ) = k
   10 CONTINUE
      IF( ( abs( a( n ) ).LE.scale1*tl ) .AND. ( in( n ).EQ.0 ) )
     $   in( n ) = n
*
      RETURN
*
*     End of DLAGTF
*
      END