1 SUBROUTINE DLAGTS( JOB, N, A, B, C, D, IN, Y, TOL, INFO )
2 *
3 * -- LAPACK auxiliary routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 INTEGER INFO, JOB, N
10 DOUBLE PRECISION TOL
11 * ..
12 * .. Array Arguments ..
13 INTEGER IN( * )
14 DOUBLE PRECISION A( * ), B( * ), C( * ), D( * ), Y( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DLAGTS may be used to solve one of the systems of equations
21 *
22 * (T - lambda*I)*x = y or (T - lambda*I)**T*x = y,
23 *
24 * where T is an n by n tridiagonal matrix, for x, following the
25 * factorization of (T - lambda*I) as
26 *
27 * (T - lambda*I) = P*L*U ,
28 *
29 * by routine DLAGTF. The choice of equation to be solved is
30 * controlled by the argument JOB, and in each case there is an option
31 * to perturb zero or very small diagonal elements of U, this option
32 * being intended for use in applications such as inverse iteration.
33 *
34 * Arguments
35 * =========
36 *
37 * JOB (input) INTEGER
38 * Specifies the job to be performed by DLAGTS as follows:
39 * = 1: The equations (T - lambda*I)x = y are to be solved,
40 * but diagonal elements of U are not to be perturbed.
41 * = -1: The equations (T - lambda*I)x = y are to be solved
42 * and, if overflow would otherwise occur, the diagonal
43 * elements of U are to be perturbed. See argument TOL
44 * below.
45 * = 2: The equations (T - lambda*I)**Tx = y are to be solved,
46 * but diagonal elements of U are not to be perturbed.
47 * = -2: The equations (T - lambda*I)**Tx = y are to be solved
48 * and, if overflow would otherwise occur, the diagonal
49 * elements of U are to be perturbed. See argument TOL
50 * below.
51 *
52 * N (input) INTEGER
53 * The order of the matrix T.
54 *
55 * A (input) DOUBLE PRECISION array, dimension (N)
56 * On entry, A must contain the diagonal elements of U as
57 * returned from DLAGTF.
58 *
59 * B (input) DOUBLE PRECISION array, dimension (N-1)
60 * On entry, B must contain the first super-diagonal elements of
61 * U as returned from DLAGTF.
62 *
63 * C (input) DOUBLE PRECISION array, dimension (N-1)
64 * On entry, C must contain the sub-diagonal elements of L as
65 * returned from DLAGTF.
66 *
67 * D (input) DOUBLE PRECISION array, dimension (N-2)
68 * On entry, D must contain the second super-diagonal elements
69 * of U as returned from DLAGTF.
70 *
71 * IN (input) INTEGER array, dimension (N)
72 * On entry, IN must contain details of the matrix P as returned
73 * from DLAGTF.
74 *
75 * Y (input/output) DOUBLE PRECISION array, dimension (N)
76 * On entry, the right hand side vector y.
77 * On exit, Y is overwritten by the solution vector x.
78 *
79 * TOL (input/output) DOUBLE PRECISION
80 * On entry, with JOB .lt. 0, TOL should be the minimum
81 * perturbation to be made to very small diagonal elements of U.
82 * TOL should normally be chosen as about eps*norm(U), where eps
83 * is the relative machine precision, but if TOL is supplied as
84 * non-positive, then it is reset to eps*max( abs( u(i,j) ) ).
85 * If JOB .gt. 0 then TOL is not referenced.
86 *
87 * On exit, TOL is changed as described above, only if TOL is
88 * non-positive on entry. Otherwise TOL is unchanged.
89 *
90 * INFO (output) INTEGER
91 * = 0 : successful exit
92 * .lt. 0: if INFO = -i, the i-th argument had an illegal value
93 * .gt. 0: overflow would occur when computing the INFO(th)
94 * element of the solution vector x. This can only occur
95 * when JOB is supplied as positive and either means
96 * that a diagonal element of U is very small, or that
97 * the elements of the right-hand side vector y are very
98 * large.
99 *
100 * =====================================================================
101 *
102 * .. Parameters ..
103 DOUBLE PRECISION ONE, ZERO
104 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
105 * ..
106 * .. Local Scalars ..
107 INTEGER K
108 DOUBLE PRECISION ABSAK, AK, BIGNUM, EPS, PERT, SFMIN, TEMP
109 * ..
110 * .. Intrinsic Functions ..
111 INTRINSIC ABS, MAX, SIGN
112 * ..
113 * .. External Functions ..
114 DOUBLE PRECISION DLAMCH
115 EXTERNAL DLAMCH
116 * ..
117 * .. External Subroutines ..
118 EXTERNAL XERBLA
119 * ..
120 * .. Executable Statements ..
121 *
122 INFO = 0
123 IF( ( ABS( JOB ).GT.2 ) .OR. ( JOB.EQ.0 ) ) THEN
124 INFO = -1
125 ELSE IF( N.LT.0 ) THEN
126 INFO = -2
127 END IF
128 IF( INFO.NE.0 ) THEN
129 CALL XERBLA( 'DLAGTS', -INFO )
130 RETURN
131 END IF
132 *
133 IF( N.EQ.0 )
134 $ RETURN
135 *
136 EPS = DLAMCH( 'Epsilon' )
137 SFMIN = DLAMCH( 'Safe minimum' )
138 BIGNUM = ONE / SFMIN
139 *
140 IF( JOB.LT.0 ) THEN
141 IF( TOL.LE.ZERO ) THEN
142 TOL = ABS( A( 1 ) )
143 IF( N.GT.1 )
144 $ TOL = MAX( TOL, ABS( A( 2 ) ), ABS( B( 1 ) ) )
145 DO 10 K = 3, N
146 TOL = MAX( TOL, ABS( A( K ) ), ABS( B( K-1 ) ),
147 $ ABS( D( K-2 ) ) )
148 10 CONTINUE
149 TOL = TOL*EPS
150 IF( TOL.EQ.ZERO )
151 $ TOL = EPS
152 END IF
153 END IF
154 *
155 IF( ABS( JOB ).EQ.1 ) THEN
156 DO 20 K = 2, N
157 IF( IN( K-1 ).EQ.0 ) THEN
158 Y( K ) = Y( K ) - C( K-1 )*Y( K-1 )
159 ELSE
160 TEMP = Y( K-1 )
161 Y( K-1 ) = Y( K )
162 Y( K ) = TEMP - C( K-1 )*Y( K )
163 END IF
164 20 CONTINUE
165 IF( JOB.EQ.1 ) THEN
166 DO 30 K = N, 1, -1
167 IF( K.LE.N-2 ) THEN
168 TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 )
169 ELSE IF( K.EQ.N-1 ) THEN
170 TEMP = Y( K ) - B( K )*Y( K+1 )
171 ELSE
172 TEMP = Y( K )
173 END IF
174 AK = A( K )
175 ABSAK = ABS( AK )
176 IF( ABSAK.LT.ONE ) THEN
177 IF( ABSAK.LT.SFMIN ) THEN
178 IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
179 $ THEN
180 INFO = K
181 RETURN
182 ELSE
183 TEMP = TEMP*BIGNUM
184 AK = AK*BIGNUM
185 END IF
186 ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
187 INFO = K
188 RETURN
189 END IF
190 END IF
191 Y( K ) = TEMP / AK
192 30 CONTINUE
193 ELSE
194 DO 50 K = N, 1, -1
195 IF( K.LE.N-2 ) THEN
196 TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 )
197 ELSE IF( K.EQ.N-1 ) THEN
198 TEMP = Y( K ) - B( K )*Y( K+1 )
199 ELSE
200 TEMP = Y( K )
201 END IF
202 AK = A( K )
203 PERT = SIGN( TOL, AK )
204 40 CONTINUE
205 ABSAK = ABS( AK )
206 IF( ABSAK.LT.ONE ) THEN
207 IF( ABSAK.LT.SFMIN ) THEN
208 IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
209 $ THEN
210 AK = AK + PERT
211 PERT = 2*PERT
212 GO TO 40
213 ELSE
214 TEMP = TEMP*BIGNUM
215 AK = AK*BIGNUM
216 END IF
217 ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
218 AK = AK + PERT
219 PERT = 2*PERT
220 GO TO 40
221 END IF
222 END IF
223 Y( K ) = TEMP / AK
224 50 CONTINUE
225 END IF
226 ELSE
227 *
228 * Come to here if JOB = 2 or -2
229 *
230 IF( JOB.EQ.2 ) THEN
231 DO 60 K = 1, N
232 IF( K.GE.3 ) THEN
233 TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 )
234 ELSE IF( K.EQ.2 ) THEN
235 TEMP = Y( K ) - B( K-1 )*Y( K-1 )
236 ELSE
237 TEMP = Y( K )
238 END IF
239 AK = A( K )
240 ABSAK = ABS( AK )
241 IF( ABSAK.LT.ONE ) THEN
242 IF( ABSAK.LT.SFMIN ) THEN
243 IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
244 $ THEN
245 INFO = K
246 RETURN
247 ELSE
248 TEMP = TEMP*BIGNUM
249 AK = AK*BIGNUM
250 END IF
251 ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
252 INFO = K
253 RETURN
254 END IF
255 END IF
256 Y( K ) = TEMP / AK
257 60 CONTINUE
258 ELSE
259 DO 80 K = 1, N
260 IF( K.GE.3 ) THEN
261 TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 )
262 ELSE IF( K.EQ.2 ) THEN
263 TEMP = Y( K ) - B( K-1 )*Y( K-1 )
264 ELSE
265 TEMP = Y( K )
266 END IF
267 AK = A( K )
268 PERT = SIGN( TOL, AK )
269 70 CONTINUE
270 ABSAK = ABS( AK )
271 IF( ABSAK.LT.ONE ) THEN
272 IF( ABSAK.LT.SFMIN ) THEN
273 IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
274 $ THEN
275 AK = AK + PERT
276 PERT = 2*PERT
277 GO TO 70
278 ELSE
279 TEMP = TEMP*BIGNUM
280 AK = AK*BIGNUM
281 END IF
282 ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
283 AK = AK + PERT
284 PERT = 2*PERT
285 GO TO 70
286 END IF
287 END IF
288 Y( K ) = TEMP / AK
289 80 CONTINUE
290 END IF
291 *
292 DO 90 K = N, 2, -1
293 IF( IN( K-1 ).EQ.0 ) THEN
294 Y( K-1 ) = Y( K-1 ) - C( K-1 )*Y( K )
295 ELSE
296 TEMP = Y( K-1 )
297 Y( K-1 ) = Y( K )
298 Y( K ) = TEMP - C( K-1 )*Y( K )
299 END IF
300 90 CONTINUE
301 END IF
302 *
303 * End of DLAGTS
304 *
305 END
2 *
3 * -- LAPACK auxiliary routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 INTEGER INFO, JOB, N
10 DOUBLE PRECISION TOL
11 * ..
12 * .. Array Arguments ..
13 INTEGER IN( * )
14 DOUBLE PRECISION A( * ), B( * ), C( * ), D( * ), Y( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DLAGTS may be used to solve one of the systems of equations
21 *
22 * (T - lambda*I)*x = y or (T - lambda*I)**T*x = y,
23 *
24 * where T is an n by n tridiagonal matrix, for x, following the
25 * factorization of (T - lambda*I) as
26 *
27 * (T - lambda*I) = P*L*U ,
28 *
29 * by routine DLAGTF. The choice of equation to be solved is
30 * controlled by the argument JOB, and in each case there is an option
31 * to perturb zero or very small diagonal elements of U, this option
32 * being intended for use in applications such as inverse iteration.
33 *
34 * Arguments
35 * =========
36 *
37 * JOB (input) INTEGER
38 * Specifies the job to be performed by DLAGTS as follows:
39 * = 1: The equations (T - lambda*I)x = y are to be solved,
40 * but diagonal elements of U are not to be perturbed.
41 * = -1: The equations (T - lambda*I)x = y are to be solved
42 * and, if overflow would otherwise occur, the diagonal
43 * elements of U are to be perturbed. See argument TOL
44 * below.
45 * = 2: The equations (T - lambda*I)**Tx = y are to be solved,
46 * but diagonal elements of U are not to be perturbed.
47 * = -2: The equations (T - lambda*I)**Tx = y are to be solved
48 * and, if overflow would otherwise occur, the diagonal
49 * elements of U are to be perturbed. See argument TOL
50 * below.
51 *
52 * N (input) INTEGER
53 * The order of the matrix T.
54 *
55 * A (input) DOUBLE PRECISION array, dimension (N)
56 * On entry, A must contain the diagonal elements of U as
57 * returned from DLAGTF.
58 *
59 * B (input) DOUBLE PRECISION array, dimension (N-1)
60 * On entry, B must contain the first super-diagonal elements of
61 * U as returned from DLAGTF.
62 *
63 * C (input) DOUBLE PRECISION array, dimension (N-1)
64 * On entry, C must contain the sub-diagonal elements of L as
65 * returned from DLAGTF.
66 *
67 * D (input) DOUBLE PRECISION array, dimension (N-2)
68 * On entry, D must contain the second super-diagonal elements
69 * of U as returned from DLAGTF.
70 *
71 * IN (input) INTEGER array, dimension (N)
72 * On entry, IN must contain details of the matrix P as returned
73 * from DLAGTF.
74 *
75 * Y (input/output) DOUBLE PRECISION array, dimension (N)
76 * On entry, the right hand side vector y.
77 * On exit, Y is overwritten by the solution vector x.
78 *
79 * TOL (input/output) DOUBLE PRECISION
80 * On entry, with JOB .lt. 0, TOL should be the minimum
81 * perturbation to be made to very small diagonal elements of U.
82 * TOL should normally be chosen as about eps*norm(U), where eps
83 * is the relative machine precision, but if TOL is supplied as
84 * non-positive, then it is reset to eps*max( abs( u(i,j) ) ).
85 * If JOB .gt. 0 then TOL is not referenced.
86 *
87 * On exit, TOL is changed as described above, only if TOL is
88 * non-positive on entry. Otherwise TOL is unchanged.
89 *
90 * INFO (output) INTEGER
91 * = 0 : successful exit
92 * .lt. 0: if INFO = -i, the i-th argument had an illegal value
93 * .gt. 0: overflow would occur when computing the INFO(th)
94 * element of the solution vector x. This can only occur
95 * when JOB is supplied as positive and either means
96 * that a diagonal element of U is very small, or that
97 * the elements of the right-hand side vector y are very
98 * large.
99 *
100 * =====================================================================
101 *
102 * .. Parameters ..
103 DOUBLE PRECISION ONE, ZERO
104 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
105 * ..
106 * .. Local Scalars ..
107 INTEGER K
108 DOUBLE PRECISION ABSAK, AK, BIGNUM, EPS, PERT, SFMIN, TEMP
109 * ..
110 * .. Intrinsic Functions ..
111 INTRINSIC ABS, MAX, SIGN
112 * ..
113 * .. External Functions ..
114 DOUBLE PRECISION DLAMCH
115 EXTERNAL DLAMCH
116 * ..
117 * .. External Subroutines ..
118 EXTERNAL XERBLA
119 * ..
120 * .. Executable Statements ..
121 *
122 INFO = 0
123 IF( ( ABS( JOB ).GT.2 ) .OR. ( JOB.EQ.0 ) ) THEN
124 INFO = -1
125 ELSE IF( N.LT.0 ) THEN
126 INFO = -2
127 END IF
128 IF( INFO.NE.0 ) THEN
129 CALL XERBLA( 'DLAGTS', -INFO )
130 RETURN
131 END IF
132 *
133 IF( N.EQ.0 )
134 $ RETURN
135 *
136 EPS = DLAMCH( 'Epsilon' )
137 SFMIN = DLAMCH( 'Safe minimum' )
138 BIGNUM = ONE / SFMIN
139 *
140 IF( JOB.LT.0 ) THEN
141 IF( TOL.LE.ZERO ) THEN
142 TOL = ABS( A( 1 ) )
143 IF( N.GT.1 )
144 $ TOL = MAX( TOL, ABS( A( 2 ) ), ABS( B( 1 ) ) )
145 DO 10 K = 3, N
146 TOL = MAX( TOL, ABS( A( K ) ), ABS( B( K-1 ) ),
147 $ ABS( D( K-2 ) ) )
148 10 CONTINUE
149 TOL = TOL*EPS
150 IF( TOL.EQ.ZERO )
151 $ TOL = EPS
152 END IF
153 END IF
154 *
155 IF( ABS( JOB ).EQ.1 ) THEN
156 DO 20 K = 2, N
157 IF( IN( K-1 ).EQ.0 ) THEN
158 Y( K ) = Y( K ) - C( K-1 )*Y( K-1 )
159 ELSE
160 TEMP = Y( K-1 )
161 Y( K-1 ) = Y( K )
162 Y( K ) = TEMP - C( K-1 )*Y( K )
163 END IF
164 20 CONTINUE
165 IF( JOB.EQ.1 ) THEN
166 DO 30 K = N, 1, -1
167 IF( K.LE.N-2 ) THEN
168 TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 )
169 ELSE IF( K.EQ.N-1 ) THEN
170 TEMP = Y( K ) - B( K )*Y( K+1 )
171 ELSE
172 TEMP = Y( K )
173 END IF
174 AK = A( K )
175 ABSAK = ABS( AK )
176 IF( ABSAK.LT.ONE ) THEN
177 IF( ABSAK.LT.SFMIN ) THEN
178 IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
179 $ THEN
180 INFO = K
181 RETURN
182 ELSE
183 TEMP = TEMP*BIGNUM
184 AK = AK*BIGNUM
185 END IF
186 ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
187 INFO = K
188 RETURN
189 END IF
190 END IF
191 Y( K ) = TEMP / AK
192 30 CONTINUE
193 ELSE
194 DO 50 K = N, 1, -1
195 IF( K.LE.N-2 ) THEN
196 TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 )
197 ELSE IF( K.EQ.N-1 ) THEN
198 TEMP = Y( K ) - B( K )*Y( K+1 )
199 ELSE
200 TEMP = Y( K )
201 END IF
202 AK = A( K )
203 PERT = SIGN( TOL, AK )
204 40 CONTINUE
205 ABSAK = ABS( AK )
206 IF( ABSAK.LT.ONE ) THEN
207 IF( ABSAK.LT.SFMIN ) THEN
208 IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
209 $ THEN
210 AK = AK + PERT
211 PERT = 2*PERT
212 GO TO 40
213 ELSE
214 TEMP = TEMP*BIGNUM
215 AK = AK*BIGNUM
216 END IF
217 ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
218 AK = AK + PERT
219 PERT = 2*PERT
220 GO TO 40
221 END IF
222 END IF
223 Y( K ) = TEMP / AK
224 50 CONTINUE
225 END IF
226 ELSE
227 *
228 * Come to here if JOB = 2 or -2
229 *
230 IF( JOB.EQ.2 ) THEN
231 DO 60 K = 1, N
232 IF( K.GE.3 ) THEN
233 TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 )
234 ELSE IF( K.EQ.2 ) THEN
235 TEMP = Y( K ) - B( K-1 )*Y( K-1 )
236 ELSE
237 TEMP = Y( K )
238 END IF
239 AK = A( K )
240 ABSAK = ABS( AK )
241 IF( ABSAK.LT.ONE ) THEN
242 IF( ABSAK.LT.SFMIN ) THEN
243 IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
244 $ THEN
245 INFO = K
246 RETURN
247 ELSE
248 TEMP = TEMP*BIGNUM
249 AK = AK*BIGNUM
250 END IF
251 ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
252 INFO = K
253 RETURN
254 END IF
255 END IF
256 Y( K ) = TEMP / AK
257 60 CONTINUE
258 ELSE
259 DO 80 K = 1, N
260 IF( K.GE.3 ) THEN
261 TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 )
262 ELSE IF( K.EQ.2 ) THEN
263 TEMP = Y( K ) - B( K-1 )*Y( K-1 )
264 ELSE
265 TEMP = Y( K )
266 END IF
267 AK = A( K )
268 PERT = SIGN( TOL, AK )
269 70 CONTINUE
270 ABSAK = ABS( AK )
271 IF( ABSAK.LT.ONE ) THEN
272 IF( ABSAK.LT.SFMIN ) THEN
273 IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
274 $ THEN
275 AK = AK + PERT
276 PERT = 2*PERT
277 GO TO 70
278 ELSE
279 TEMP = TEMP*BIGNUM
280 AK = AK*BIGNUM
281 END IF
282 ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
283 AK = AK + PERT
284 PERT = 2*PERT
285 GO TO 70
286 END IF
287 END IF
288 Y( K ) = TEMP / AK
289 80 CONTINUE
290 END IF
291 *
292 DO 90 K = N, 2, -1
293 IF( IN( K-1 ).EQ.0 ) THEN
294 Y( K-1 ) = Y( K-1 ) - C( K-1 )*Y( K )
295 ELSE
296 TEMP = Y( K-1 )
297 Y( K-1 ) = Y( K )
298 Y( K ) = TEMP - C( K-1 )*Y( K )
299 END IF
300 90 CONTINUE
301 END IF
302 *
303 * End of DLAGTS
304 *
305 END