001:       SUBROUTINE SLAGTF( N, A, LAMBDA, B, C, TOL, D, IN, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
005: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       INTEGER            INFO, N
010:       REAL               LAMBDA, TOL
011: *     ..
012: *     .. Array Arguments ..
013:       INTEGER            IN( * )
014:       REAL               A( * ), B( * ), C( * ), D( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  SLAGTF factorizes the matrix (T - lambda*I), where T is an n by n
021: *  tridiagonal matrix and lambda is a scalar, as
022: *
023: *     T - lambda*I = PLU,
024: *
025: *  where P is a permutation matrix, L is a unit lower tridiagonal matrix
026: *  with at most one non-zero sub-diagonal elements per column and U is
027: *  an upper triangular matrix with at most two non-zero super-diagonal
028: *  elements per column.
029: *
030: *  The factorization is obtained by Gaussian elimination with partial
031: *  pivoting and implicit row scaling.
032: *
033: *  The parameter LAMBDA is included in the routine so that SLAGTF may
034: *  be used, in conjunction with SLAGTS, to obtain eigenvectors of T by
035: *  inverse iteration.
036: *
037: *  Arguments
038: *  =========
039: *
040: *  N       (input) INTEGER
041: *          The order of the matrix T.
042: *
043: *  A       (input/output) REAL array, dimension (N)
044: *          On entry, A must contain the diagonal elements of T.
045: *
046: *          On exit, A is overwritten by the n diagonal elements of the
047: *          upper triangular matrix U of the factorization of T.
048: *
049: *  LAMBDA  (input) REAL
050: *          On entry, the scalar lambda.
051: *
052: *  B       (input/output) REAL array, dimension (N-1)
053: *          On entry, B must contain the (n-1) super-diagonal elements of
054: *          T.
055: *
056: *          On exit, B is overwritten by the (n-1) super-diagonal
057: *          elements of the matrix U of the factorization of T.
058: *
059: *  C       (input/output) REAL array, dimension (N-1)
060: *          On entry, C must contain the (n-1) sub-diagonal elements of
061: *          T.
062: *
063: *          On exit, C is overwritten by the (n-1) sub-diagonal elements
064: *          of the matrix L of the factorization of T.
065: *
066: *  TOL     (input) REAL
067: *          On entry, a relative tolerance used to indicate whether or
068: *          not the matrix (T - lambda*I) is nearly singular. TOL should
069: *          normally be chose as approximately the largest relative error
070: *          in the elements of T. For example, if the elements of T are
071: *          correct to about 4 significant figures, then TOL should be
072: *          set to about 5*10**(-4). If TOL is supplied as less than eps,
073: *          where eps is the relative machine precision, then the value
074: *          eps is used in place of TOL.
075: *
076: *  D       (output) REAL array, dimension (N-2)
077: *          On exit, D is overwritten by the (n-2) second super-diagonal
078: *          elements of the matrix U of the factorization of T.
079: *
080: *  IN      (output) INTEGER array, dimension (N)
081: *          On exit, IN contains details of the permutation matrix P. If
082: *          an interchange occurred at the kth step of the elimination,
083: *          then IN(k) = 1, otherwise IN(k) = 0. The element IN(n)
084: *          returns the smallest positive integer j such that
085: *
086: *             abs( u(j,j) ).le. norm( (T - lambda*I)(j) )*TOL,
087: *
088: *          where norm( A(j) ) denotes the sum of the absolute values of
089: *          the jth row of the matrix A. If no such j exists then IN(n)
090: *          is returned as zero. If IN(n) is returned as positive, then a
091: *          diagonal element of U is small, indicating that
092: *          (T - lambda*I) is singular or nearly singular,
093: *
094: *  INFO    (output) INTEGER
095: *          = 0   : successful exit
096: *          .lt. 0: if INFO = -k, the kth argument had an illegal value
097: *
098: * =====================================================================
099: *
100: *     .. Parameters ..
101:       REAL               ZERO
102:       PARAMETER          ( ZERO = 0.0E+0 )
103: *     ..
104: *     .. Local Scalars ..
105:       INTEGER            K
106:       REAL               EPS, MULT, PIV1, PIV2, SCALE1, SCALE2, TEMP, TL
107: *     ..
108: *     .. Intrinsic Functions ..
109:       INTRINSIC          ABS, MAX
110: *     ..
111: *     .. External Functions ..
112:       REAL               SLAMCH
113:       EXTERNAL           SLAMCH
114: *     ..
115: *     .. External Subroutines ..
116:       EXTERNAL           XERBLA
117: *     ..
118: *     .. Executable Statements ..
119: *
120:       INFO = 0
121:       IF( N.LT.0 ) THEN
122:          INFO = -1
123:          CALL XERBLA( 'SLAGTF', -INFO )
124:          RETURN
125:       END IF
126: *
127:       IF( N.EQ.0 )
128:      $   RETURN
129: *
130:       A( 1 ) = A( 1 ) - LAMBDA
131:       IN( N ) = 0
132:       IF( N.EQ.1 ) THEN
133:          IF( A( 1 ).EQ.ZERO )
134:      $      IN( 1 ) = 1
135:          RETURN
136:       END IF
137: *
138:       EPS = SLAMCH( 'Epsilon' )
139: *
140:       TL = MAX( TOL, EPS )
141:       SCALE1 = ABS( A( 1 ) ) + ABS( B( 1 ) )
142:       DO 10 K = 1, N - 1
143:          A( K+1 ) = A( K+1 ) - LAMBDA
144:          SCALE2 = ABS( C( K ) ) + ABS( A( K+1 ) )
145:          IF( K.LT.( N-1 ) )
146:      $      SCALE2 = SCALE2 + ABS( B( K+1 ) )
147:          IF( A( K ).EQ.ZERO ) THEN
148:             PIV1 = ZERO
149:          ELSE
150:             PIV1 = ABS( A( K ) ) / SCALE1
151:          END IF
152:          IF( C( K ).EQ.ZERO ) THEN
153:             IN( K ) = 0
154:             PIV2 = ZERO
155:             SCALE1 = SCALE2
156:             IF( K.LT.( N-1 ) )
157:      $         D( K ) = ZERO
158:          ELSE
159:             PIV2 = ABS( C( K ) ) / SCALE2
160:             IF( PIV2.LE.PIV1 ) THEN
161:                IN( K ) = 0
162:                SCALE1 = SCALE2
163:                C( K ) = C( K ) / A( K )
164:                A( K+1 ) = A( K+1 ) - C( K )*B( K )
165:                IF( K.LT.( N-1 ) )
166:      $            D( K ) = ZERO
167:             ELSE
168:                IN( K ) = 1
169:                MULT = A( K ) / C( K )
170:                A( K ) = C( K )
171:                TEMP = A( K+1 )
172:                A( K+1 ) = B( K ) - MULT*TEMP
173:                IF( K.LT.( N-1 ) ) THEN
174:                   D( K ) = B( K+1 )
175:                   B( K+1 ) = -MULT*D( K )
176:                END IF
177:                B( K ) = TEMP
178:                C( K ) = MULT
179:             END IF
180:          END IF
181:          IF( ( MAX( PIV1, PIV2 ).LE.TL ) .AND. ( IN( N ).EQ.0 ) )
182:      $      IN( N ) = K
183:    10 CONTINUE
184:       IF( ( ABS( A( N ) ).LE.SCALE1*TL ) .AND. ( IN( N ).EQ.0 ) )
185:      $   IN( N ) = N
186: *
187:       RETURN
188: *
189: *     End of SLAGTF
190: *
191:       END
192: