LAPACK 3.3.1 Linear Algebra PACKage

# dlagtf.f

Go to the documentation of this file.
```00001       SUBROUTINE DLAGTF( N, A, LAMBDA, B, C, TOL, D, IN, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            INFO, N
00010       DOUBLE PRECISION   LAMBDA, TOL
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            IN( * )
00014       DOUBLE PRECISION   A( * ), B( * ), C( * ), D( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  DLAGTF factorizes the matrix (T - lambda*I), where T is an n by n
00021 *  tridiagonal matrix and lambda is a scalar, as
00022 *
00023 *     T - lambda*I = PLU,
00024 *
00025 *  where P is a permutation matrix, L is a unit lower tridiagonal matrix
00026 *  with at most one non-zero sub-diagonal elements per column and U is
00027 *  an upper triangular matrix with at most two non-zero super-diagonal
00028 *  elements per column.
00029 *
00030 *  The factorization is obtained by Gaussian elimination with partial
00031 *  pivoting and implicit row scaling.
00032 *
00033 *  The parameter LAMBDA is included in the routine so that DLAGTF may
00034 *  be used, in conjunction with DLAGTS, to obtain eigenvectors of T by
00035 *  inverse iteration.
00036 *
00037 *  Arguments
00038 *  =========
00039 *
00040 *  N       (input) INTEGER
00041 *          The order of the matrix T.
00042 *
00043 *  A       (input/output) DOUBLE PRECISION array, dimension (N)
00044 *          On entry, A must contain the diagonal elements of T.
00045 *
00046 *          On exit, A is overwritten by the n diagonal elements of the
00047 *          upper triangular matrix U of the factorization of T.
00048 *
00049 *  LAMBDA  (input) DOUBLE PRECISION
00050 *          On entry, the scalar lambda.
00051 *
00052 *  B       (input/output) DOUBLE PRECISION array, dimension (N-1)
00053 *          On entry, B must contain the (n-1) super-diagonal elements of
00054 *          T.
00055 *
00056 *          On exit, B is overwritten by the (n-1) super-diagonal
00057 *          elements of the matrix U of the factorization of T.
00058 *
00059 *  C       (input/output) DOUBLE PRECISION array, dimension (N-1)
00060 *          On entry, C must contain the (n-1) sub-diagonal elements of
00061 *          T.
00062 *
00063 *          On exit, C is overwritten by the (n-1) sub-diagonal elements
00064 *          of the matrix L of the factorization of T.
00065 *
00066 *  TOL     (input) DOUBLE PRECISION
00067 *          On entry, a relative tolerance used to indicate whether or
00068 *          not the matrix (T - lambda*I) is nearly singular. TOL should
00069 *          normally be chose as approximately the largest relative error
00070 *          in the elements of T. For example, if the elements of T are
00071 *          correct to about 4 significant figures, then TOL should be
00072 *          set to about 5*10**(-4). If TOL is supplied as less than eps,
00073 *          where eps is the relative machine precision, then the value
00074 *          eps is used in place of TOL.
00075 *
00076 *  D       (output) DOUBLE PRECISION array, dimension (N-2)
00077 *          On exit, D is overwritten by the (n-2) second super-diagonal
00078 *          elements of the matrix U of the factorization of T.
00079 *
00080 *  IN      (output) INTEGER array, dimension (N)
00081 *          On exit, IN contains details of the permutation matrix P. If
00082 *          an interchange occurred at the kth step of the elimination,
00083 *          then IN(k) = 1, otherwise IN(k) = 0. The element IN(n)
00084 *          returns the smallest positive integer j such that
00085 *
00086 *             abs( u(j,j) ).le. norm( (T - lambda*I)(j) )*TOL,
00087 *
00088 *          where norm( A(j) ) denotes the sum of the absolute values of
00089 *          the jth row of the matrix A. If no such j exists then IN(n)
00090 *          is returned as zero. If IN(n) is returned as positive, then a
00091 *          diagonal element of U is small, indicating that
00092 *          (T - lambda*I) is singular or nearly singular,
00093 *
00094 *  INFO    (output) INTEGER
00095 *          = 0   : successful exit
00096 *          .lt. 0: if INFO = -k, the kth argument had an illegal value
00097 *
00098 * =====================================================================
00099 *
00100 *     .. Parameters ..
00101       DOUBLE PRECISION   ZERO
00102       PARAMETER          ( ZERO = 0.0D+0 )
00103 *     ..
00104 *     .. Local Scalars ..
00105       INTEGER            K
00106       DOUBLE PRECISION   EPS, MULT, PIV1, PIV2, SCALE1, SCALE2, TEMP, TL
00107 *     ..
00108 *     .. Intrinsic Functions ..
00109       INTRINSIC          ABS, MAX
00110 *     ..
00111 *     .. External Functions ..
00112       DOUBLE PRECISION   DLAMCH
00113       EXTERNAL           DLAMCH
00114 *     ..
00115 *     .. External Subroutines ..
00116       EXTERNAL           XERBLA
00117 *     ..
00118 *     .. Executable Statements ..
00119 *
00120       INFO = 0
00121       IF( N.LT.0 ) THEN
00122          INFO = -1
00123          CALL XERBLA( 'DLAGTF', -INFO )
00124          RETURN
00125       END IF
00126 *
00127       IF( N.EQ.0 )
00128      \$   RETURN
00129 *
00130       A( 1 ) = A( 1 ) - LAMBDA
00131       IN( N ) = 0
00132       IF( N.EQ.1 ) THEN
00133          IF( A( 1 ).EQ.ZERO )
00134      \$      IN( 1 ) = 1
00135          RETURN
00136       END IF
00137 *
00138       EPS = DLAMCH( 'Epsilon' )
00139 *
00140       TL = MAX( TOL, EPS )
00141       SCALE1 = ABS( A( 1 ) ) + ABS( B( 1 ) )
00142       DO 10 K = 1, N - 1
00143          A( K+1 ) = A( K+1 ) - LAMBDA
00144          SCALE2 = ABS( C( K ) ) + ABS( A( K+1 ) )
00145          IF( K.LT.( N-1 ) )
00146      \$      SCALE2 = SCALE2 + ABS( B( K+1 ) )
00147          IF( A( K ).EQ.ZERO ) THEN
00148             PIV1 = ZERO
00149          ELSE
00150             PIV1 = ABS( A( K ) ) / SCALE1
00151          END IF
00152          IF( C( K ).EQ.ZERO ) THEN
00153             IN( K ) = 0
00154             PIV2 = ZERO
00155             SCALE1 = SCALE2
00156             IF( K.LT.( N-1 ) )
00157      \$         D( K ) = ZERO
00158          ELSE
00159             PIV2 = ABS( C( K ) ) / SCALE2
00160             IF( PIV2.LE.PIV1 ) THEN
00161                IN( K ) = 0
00162                SCALE1 = SCALE2
00163                C( K ) = C( K ) / A( K )
00164                A( K+1 ) = A( K+1 ) - C( K )*B( K )
00165                IF( K.LT.( N-1 ) )
00166      \$            D( K ) = ZERO
00167             ELSE
00168                IN( K ) = 1
00169                MULT = A( K ) / C( K )
00170                A( K ) = C( K )
00171                TEMP = A( K+1 )
00172                A( K+1 ) = B( K ) - MULT*TEMP
00173                IF( K.LT.( N-1 ) ) THEN
00174                   D( K ) = B( K+1 )
00175                   B( K+1 ) = -MULT*D( K )
00176                END IF
00177                B( K ) = TEMP
00178                C( K ) = MULT
00179             END IF
00180          END IF
00181          IF( ( MAX( PIV1, PIV2 ).LE.TL ) .AND. ( IN( N ).EQ.0 ) )
00182      \$      IN( N ) = K
00183    10 CONTINUE
00184       IF( ( ABS( A( N ) ).LE.SCALE1*TL ) .AND. ( IN( N ).EQ.0 ) )
00185      \$   IN( N ) = N
00186 *
00187       RETURN
00188 *
00189 *     End of DLAGTF
00190 *
00191       END
```