1 SUBROUTINE DLAGTF( N, A, LAMBDA, B, C, TOL, D, IN, INFO )
2 *
3 * -- LAPACK routine (version 3.2) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER INFO, N
10 DOUBLE PRECISION LAMBDA, TOL
11 * ..
12 * .. Array Arguments ..
13 INTEGER IN( * )
14 DOUBLE PRECISION A( * ), B( * ), C( * ), D( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DLAGTF factorizes the matrix (T - lambda*I), where T is an n by n
21 * tridiagonal matrix and lambda is a scalar, as
22 *
23 * T - lambda*I = PLU,
24 *
25 * where P is a permutation matrix, L is a unit lower tridiagonal matrix
26 * with at most one non-zero sub-diagonal elements per column and U is
27 * an upper triangular matrix with at most two non-zero super-diagonal
28 * elements per column.
29 *
30 * The factorization is obtained by Gaussian elimination with partial
31 * pivoting and implicit row scaling.
32 *
33 * The parameter LAMBDA is included in the routine so that DLAGTF may
34 * be used, in conjunction with DLAGTS, to obtain eigenvectors of T by
35 * inverse iteration.
36 *
37 * Arguments
38 * =========
39 *
40 * N (input) INTEGER
41 * The order of the matrix T.
42 *
43 * A (input/output) DOUBLE PRECISION array, dimension (N)
44 * On entry, A must contain the diagonal elements of T.
45 *
46 * On exit, A is overwritten by the n diagonal elements of the
47 * upper triangular matrix U of the factorization of T.
48 *
49 * LAMBDA (input) DOUBLE PRECISION
50 * On entry, the scalar lambda.
51 *
52 * B (input/output) DOUBLE PRECISION array, dimension (N-1)
53 * On entry, B must contain the (n-1) super-diagonal elements of
54 * T.
55 *
56 * On exit, B is overwritten by the (n-1) super-diagonal
57 * elements of the matrix U of the factorization of T.
58 *
59 * C (input/output) DOUBLE PRECISION array, dimension (N-1)
60 * On entry, C must contain the (n-1) sub-diagonal elements of
61 * T.
62 *
63 * On exit, C is overwritten by the (n-1) sub-diagonal elements
64 * of the matrix L of the factorization of T.
65 *
66 * TOL (input) DOUBLE PRECISION
67 * On entry, a relative tolerance used to indicate whether or
68 * not the matrix (T - lambda*I) is nearly singular. TOL should
69 * normally be chose as approximately the largest relative error
70 * in the elements of T. For example, if the elements of T are
71 * correct to about 4 significant figures, then TOL should be
72 * set to about 5*10**(-4). If TOL is supplied as less than eps,
73 * where eps is the relative machine precision, then the value
74 * eps is used in place of TOL.
75 *
76 * D (output) DOUBLE PRECISION array, dimension (N-2)
77 * On exit, D is overwritten by the (n-2) second super-diagonal
78 * elements of the matrix U of the factorization of T.
79 *
80 * IN (output) INTEGER array, dimension (N)
81 * On exit, IN contains details of the permutation matrix P. If
82 * an interchange occurred at the kth step of the elimination,
83 * then IN(k) = 1, otherwise IN(k) = 0. The element IN(n)
84 * returns the smallest positive integer j such that
85 *
86 * abs( u(j,j) ).le. norm( (T - lambda*I)(j) )*TOL,
87 *
88 * where norm( A(j) ) denotes the sum of the absolute values of
89 * the jth row of the matrix A. If no such j exists then IN(n)
90 * is returned as zero. If IN(n) is returned as positive, then a
91 * diagonal element of U is small, indicating that
92 * (T - lambda*I) is singular or nearly singular,
93 *
94 * INFO (output) INTEGER
95 * = 0 : successful exit
96 * .lt. 0: if INFO = -k, the kth argument had an illegal value
97 *
98 * =====================================================================
99 *
100 * .. Parameters ..
101 DOUBLE PRECISION ZERO
102 PARAMETER ( ZERO = 0.0D+0 )
103 * ..
104 * .. Local Scalars ..
105 INTEGER K
106 DOUBLE PRECISION EPS, MULT, PIV1, PIV2, SCALE1, SCALE2, TEMP, TL
107 * ..
108 * .. Intrinsic Functions ..
109 INTRINSIC ABS, MAX
110 * ..
111 * .. External Functions ..
112 DOUBLE PRECISION DLAMCH
113 EXTERNAL DLAMCH
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( 'DLAGTF', -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 = DLAMCH( '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 DLAGTF
190 *
191 END
2 *
3 * -- LAPACK routine (version 3.2) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER INFO, N
10 DOUBLE PRECISION LAMBDA, TOL
11 * ..
12 * .. Array Arguments ..
13 INTEGER IN( * )
14 DOUBLE PRECISION A( * ), B( * ), C( * ), D( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DLAGTF factorizes the matrix (T - lambda*I), where T is an n by n
21 * tridiagonal matrix and lambda is a scalar, as
22 *
23 * T - lambda*I = PLU,
24 *
25 * where P is a permutation matrix, L is a unit lower tridiagonal matrix
26 * with at most one non-zero sub-diagonal elements per column and U is
27 * an upper triangular matrix with at most two non-zero super-diagonal
28 * elements per column.
29 *
30 * The factorization is obtained by Gaussian elimination with partial
31 * pivoting and implicit row scaling.
32 *
33 * The parameter LAMBDA is included in the routine so that DLAGTF may
34 * be used, in conjunction with DLAGTS, to obtain eigenvectors of T by
35 * inverse iteration.
36 *
37 * Arguments
38 * =========
39 *
40 * N (input) INTEGER
41 * The order of the matrix T.
42 *
43 * A (input/output) DOUBLE PRECISION array, dimension (N)
44 * On entry, A must contain the diagonal elements of T.
45 *
46 * On exit, A is overwritten by the n diagonal elements of the
47 * upper triangular matrix U of the factorization of T.
48 *
49 * LAMBDA (input) DOUBLE PRECISION
50 * On entry, the scalar lambda.
51 *
52 * B (input/output) DOUBLE PRECISION array, dimension (N-1)
53 * On entry, B must contain the (n-1) super-diagonal elements of
54 * T.
55 *
56 * On exit, B is overwritten by the (n-1) super-diagonal
57 * elements of the matrix U of the factorization of T.
58 *
59 * C (input/output) DOUBLE PRECISION array, dimension (N-1)
60 * On entry, C must contain the (n-1) sub-diagonal elements of
61 * T.
62 *
63 * On exit, C is overwritten by the (n-1) sub-diagonal elements
64 * of the matrix L of the factorization of T.
65 *
66 * TOL (input) DOUBLE PRECISION
67 * On entry, a relative tolerance used to indicate whether or
68 * not the matrix (T - lambda*I) is nearly singular. TOL should
69 * normally be chose as approximately the largest relative error
70 * in the elements of T. For example, if the elements of T are
71 * correct to about 4 significant figures, then TOL should be
72 * set to about 5*10**(-4). If TOL is supplied as less than eps,
73 * where eps is the relative machine precision, then the value
74 * eps is used in place of TOL.
75 *
76 * D (output) DOUBLE PRECISION array, dimension (N-2)
77 * On exit, D is overwritten by the (n-2) second super-diagonal
78 * elements of the matrix U of the factorization of T.
79 *
80 * IN (output) INTEGER array, dimension (N)
81 * On exit, IN contains details of the permutation matrix P. If
82 * an interchange occurred at the kth step of the elimination,
83 * then IN(k) = 1, otherwise IN(k) = 0. The element IN(n)
84 * returns the smallest positive integer j such that
85 *
86 * abs( u(j,j) ).le. norm( (T - lambda*I)(j) )*TOL,
87 *
88 * where norm( A(j) ) denotes the sum of the absolute values of
89 * the jth row of the matrix A. If no such j exists then IN(n)
90 * is returned as zero. If IN(n) is returned as positive, then a
91 * diagonal element of U is small, indicating that
92 * (T - lambda*I) is singular or nearly singular,
93 *
94 * INFO (output) INTEGER
95 * = 0 : successful exit
96 * .lt. 0: if INFO = -k, the kth argument had an illegal value
97 *
98 * =====================================================================
99 *
100 * .. Parameters ..
101 DOUBLE PRECISION ZERO
102 PARAMETER ( ZERO = 0.0D+0 )
103 * ..
104 * .. Local Scalars ..
105 INTEGER K
106 DOUBLE PRECISION EPS, MULT, PIV1, PIV2, SCALE1, SCALE2, TEMP, TL
107 * ..
108 * .. Intrinsic Functions ..
109 INTRINSIC ABS, MAX
110 * ..
111 * .. External Functions ..
112 DOUBLE PRECISION DLAMCH
113 EXTERNAL DLAMCH
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( 'DLAGTF', -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 = DLAMCH( '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 DLAGTF
190 *
191 END