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
1.7.2