LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
slagts.f
Go to the documentation of this file.
1*> \brief \b SLAGTS solves the system of equations (T-λI)x = y
2*> or (T-λI)^Tx = y, where T is a general tridiagonal matrix
3*> and λ a scalar, using the LU factorization computed by slagtf.
4*
5* =========== DOCUMENTATION ===========
6*
7* Online html documentation available at
8* http://www.netlib.org/lapack/explore-html/
9*
10*> Download SLAGTS + dependencies
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slagts.f">
12*> [TGZ]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slagts.f">
14*> [ZIP]</a>
15*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slagts.f">
16*> [TXT]</a>
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SLAGTS( JOB, N, A, B, C, D, IN, Y, TOL, INFO )
22*
23* .. Scalar Arguments ..
24* INTEGER INFO, JOB, N
25* REAL TOL
26* ..
27* .. Array Arguments ..
28* INTEGER IN( * )
29* REAL A( * ), B( * ), C( * ), D( * ), Y( * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> SLAGTS may be used to solve one of the systems of equations
39*>
40*> (T - lambda*I)*x = y or (T - lambda*I)**T*x = y,
41*>
42*> where T is an n by n tridiagonal matrix, for x, following the
43*> factorization of (T - lambda*I) as
44*>
45*> (T - lambda*I) = P*L*U ,
46*>
47*> by routine SLAGTF. The choice of equation to be solved is
48*> controlled by the argument JOB, and in each case there is an option
49*> to perturb zero or very small diagonal elements of U, this option
50*> being intended for use in applications such as inverse iteration.
51*> \endverbatim
52*
53* Arguments:
54* ==========
55*
56*> \param[in] JOB
57*> \verbatim
58*> JOB is INTEGER
59*> Specifies the job to be performed by SLAGTS as follows:
60*> = 1: The equations (T - lambda*I)x = y are to be solved,
61*> but diagonal elements of U are not to be perturbed.
62*> = -1: The equations (T - lambda*I)x = y are to be solved
63*> and, if overflow would otherwise occur, the diagonal
64*> elements of U are to be perturbed. See argument TOL
65*> below.
66*> = 2: The equations (T - lambda*I)**Tx = y are to be solved,
67*> but diagonal elements of U are not to be perturbed.
68*> = -2: The equations (T - lambda*I)**Tx = y are to be solved
69*> and, if overflow would otherwise occur, the diagonal
70*> elements of U are to be perturbed. See argument TOL
71*> below.
72*> \endverbatim
73*>
74*> \param[in] N
75*> \verbatim
76*> N is INTEGER
77*> The order of the matrix T.
78*> \endverbatim
79*>
80*> \param[in] A
81*> \verbatim
82*> A is REAL array, dimension (N)
83*> On entry, A must contain the diagonal elements of U as
84*> returned from SLAGTF.
85*> \endverbatim
86*>
87*> \param[in] B
88*> \verbatim
89*> B is REAL array, dimension (N-1)
90*> On entry, B must contain the first super-diagonal elements of
91*> U as returned from SLAGTF.
92*> \endverbatim
93*>
94*> \param[in] C
95*> \verbatim
96*> C is REAL array, dimension (N-1)
97*> On entry, C must contain the sub-diagonal elements of L as
98*> returned from SLAGTF.
99*> \endverbatim
100*>
101*> \param[in] D
102*> \verbatim
103*> D is REAL array, dimension (N-2)
104*> On entry, D must contain the second super-diagonal elements
105*> of U as returned from SLAGTF.
106*> \endverbatim
107*>
108*> \param[in] IN
109*> \verbatim
110*> IN is INTEGER array, dimension (N)
111*> On entry, IN must contain details of the matrix P as returned
112*> from SLAGTF.
113*> \endverbatim
114*>
115*> \param[in,out] Y
116*> \verbatim
117*> Y is REAL array, dimension (N)
118*> On entry, the right hand side vector y.
119*> On exit, Y is overwritten by the solution vector x.
120*> \endverbatim
121*>
122*> \param[in,out] TOL
123*> \verbatim
124*> TOL is REAL
125*> On entry, with JOB < 0, TOL should be the minimum
126*> perturbation to be made to very small diagonal elements of U.
127*> TOL should normally be chosen as about eps*norm(U), where eps
128*> is the relative machine precision, but if TOL is supplied as
129*> non-positive, then it is reset to eps*max( abs( u(i,j) ) ).
130*> If JOB > 0 then TOL is not referenced.
131*>
132*> On exit, TOL is changed as described above, only if TOL is
133*> non-positive on entry. Otherwise TOL is unchanged.
134*> \endverbatim
135*>
136*> \param[out] INFO
137*> \verbatim
138*> INFO is INTEGER
139*> = 0: successful exit
140*> < 0: if INFO = -i, the i-th argument had an illegal value
141*> > 0: overflow would occur when computing the INFO(th)
142*> element of the solution vector x. This can only occur
143*> when JOB is supplied as positive and either means
144*> that a diagonal element of U is very small, or that
145*> the elements of the right-hand side vector y are very
146*> large.
147*> \endverbatim
148*
149* Authors:
150* ========
151*
152*> \author Univ. of Tennessee
153*> \author Univ. of California Berkeley
154*> \author Univ. of Colorado Denver
155*> \author NAG Ltd.
156*
157*> \ingroup lagts
158*
159* =====================================================================
160 SUBROUTINE slagts( JOB, N, A, B, C, D, IN, Y, TOL, INFO )
161*
162* -- LAPACK auxiliary routine --
163* -- LAPACK is a software package provided by Univ. of Tennessee, --
164* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
165*
166* .. Scalar Arguments ..
167 INTEGER INFO, JOB, N
168 REAL TOL
169* ..
170* .. Array Arguments ..
171 INTEGER IN( * )
172 REAL A( * ), B( * ), C( * ), D( * ), Y( * )
173* ..
174*
175* =====================================================================
176*
177* .. Parameters ..
178 REAL ONE, ZERO
179 parameter( one = 1.0e+0, zero = 0.0e+0 )
180* ..
181* .. Local Scalars ..
182 INTEGER K
183 REAL ABSAK, AK, BIGNUM, EPS, PERT, SFMIN, TEMP
184* ..
185* .. Intrinsic Functions ..
186 INTRINSIC abs, max, sign
187* ..
188* .. External Functions ..
189 REAL SLAMCH
190 EXTERNAL slamch
191* ..
192* .. External Subroutines ..
193 EXTERNAL xerbla
194* ..
195* .. Executable Statements ..
196*
197 info = 0
198 IF( ( abs( job ).GT.2 ) .OR. ( job.EQ.0 ) ) THEN
199 info = -1
200 ELSE IF( n.LT.0 ) THEN
201 info = -2
202 END IF
203 IF( info.NE.0 ) THEN
204 CALL xerbla( 'SLAGTS', -info )
205 RETURN
206 END IF
207*
208 IF( n.EQ.0 )
209 $ RETURN
210*
211 eps = slamch( 'Epsilon' )
212 sfmin = slamch( 'Safe minimum' )
213 bignum = one / sfmin
214*
215 IF( job.LT.0 ) THEN
216 IF( tol.LE.zero ) THEN
217 tol = abs( a( 1 ) )
218 IF( n.GT.1 )
219 $ tol = max( tol, abs( a( 2 ) ), abs( b( 1 ) ) )
220 DO 10 k = 3, n
221 tol = max( tol, abs( a( k ) ), abs( b( k-1 ) ),
222 $ abs( d( k-2 ) ) )
223 10 CONTINUE
224 tol = tol*eps
225 IF( tol.EQ.zero )
226 $ tol = eps
227 END IF
228 END IF
229*
230 IF( abs( job ).EQ.1 ) THEN
231 DO 20 k = 2, n
232 IF( in( k-1 ).EQ.0 ) THEN
233 y( k ) = y( k ) - c( k-1 )*y( k-1 )
234 ELSE
235 temp = y( k-1 )
236 y( k-1 ) = y( k )
237 y( k ) = temp - c( k-1 )*y( k )
238 END IF
239 20 CONTINUE
240 IF( job.EQ.1 ) THEN
241 DO 30 k = n, 1, -1
242 IF( k.LE.n-2 ) THEN
243 temp = y( k ) - b( k )*y( k+1 ) - d( k )*y( k+2 )
244 ELSE IF( k.EQ.n-1 ) THEN
245 temp = y( k ) - b( k )*y( k+1 )
246 ELSE
247 temp = y( k )
248 END IF
249 ak = a( k )
250 absak = abs( ak )
251 IF( absak.LT.one ) THEN
252 IF( absak.LT.sfmin ) THEN
253 IF( absak.EQ.zero .OR. abs( temp )*sfmin.GT.absak )
254 $ THEN
255 info = k
256 RETURN
257 ELSE
258 temp = temp*bignum
259 ak = ak*bignum
260 END IF
261 ELSE IF( abs( temp ).GT.absak*bignum ) THEN
262 info = k
263 RETURN
264 END IF
265 END IF
266 y( k ) = temp / ak
267 30 CONTINUE
268 ELSE
269 DO 50 k = n, 1, -1
270 IF( k.LE.n-2 ) THEN
271 temp = y( k ) - b( k )*y( k+1 ) - d( k )*y( k+2 )
272 ELSE IF( k.EQ.n-1 ) THEN
273 temp = y( k ) - b( k )*y( k+1 )
274 ELSE
275 temp = y( k )
276 END IF
277 ak = a( k )
278 pert = sign( tol, ak )
279 40 CONTINUE
280 absak = abs( ak )
281 IF( absak.LT.one ) THEN
282 IF( absak.LT.sfmin ) THEN
283 IF( absak.EQ.zero .OR. abs( temp )*sfmin.GT.absak )
284 $ THEN
285 ak = ak + pert
286 pert = 2*pert
287 GO TO 40
288 ELSE
289 temp = temp*bignum
290 ak = ak*bignum
291 END IF
292 ELSE IF( abs( temp ).GT.absak*bignum ) THEN
293 ak = ak + pert
294 pert = 2*pert
295 GO TO 40
296 END IF
297 END IF
298 y( k ) = temp / ak
299 50 CONTINUE
300 END IF
301 ELSE
302*
303* Come to here if JOB = 2 or -2
304*
305 IF( job.EQ.2 ) THEN
306 DO 60 k = 1, n
307 IF( k.GE.3 ) THEN
308 temp = y( k ) - b( k-1 )*y( k-1 ) - d( k-2 )*y( k-2 )
309 ELSE IF( k.EQ.2 ) THEN
310 temp = y( k ) - b( k-1 )*y( k-1 )
311 ELSE
312 temp = y( k )
313 END IF
314 ak = a( k )
315 absak = abs( ak )
316 IF( absak.LT.one ) THEN
317 IF( absak.LT.sfmin ) THEN
318 IF( absak.EQ.zero .OR. abs( temp )*sfmin.GT.absak )
319 $ THEN
320 info = k
321 RETURN
322 ELSE
323 temp = temp*bignum
324 ak = ak*bignum
325 END IF
326 ELSE IF( abs( temp ).GT.absak*bignum ) THEN
327 info = k
328 RETURN
329 END IF
330 END IF
331 y( k ) = temp / ak
332 60 CONTINUE
333 ELSE
334 DO 80 k = 1, n
335 IF( k.GE.3 ) THEN
336 temp = y( k ) - b( k-1 )*y( k-1 ) - d( k-2 )*y( k-2 )
337 ELSE IF( k.EQ.2 ) THEN
338 temp = y( k ) - b( k-1 )*y( k-1 )
339 ELSE
340 temp = y( k )
341 END IF
342 ak = a( k )
343 pert = sign( tol, ak )
344 70 CONTINUE
345 absak = abs( ak )
346 IF( absak.LT.one ) THEN
347 IF( absak.LT.sfmin ) THEN
348 IF( absak.EQ.zero .OR. abs( temp )*sfmin.GT.absak )
349 $ THEN
350 ak = ak + pert
351 pert = 2*pert
352 GO TO 70
353 ELSE
354 temp = temp*bignum
355 ak = ak*bignum
356 END IF
357 ELSE IF( abs( temp ).GT.absak*bignum ) THEN
358 ak = ak + pert
359 pert = 2*pert
360 GO TO 70
361 END IF
362 END IF
363 y( k ) = temp / ak
364 80 CONTINUE
365 END IF
366*
367 DO 90 k = n, 2, -1
368 IF( in( k-1 ).EQ.0 ) THEN
369 y( k-1 ) = y( k-1 ) - c( k-1 )*y( k )
370 ELSE
371 temp = y( k-1 )
372 y( k-1 ) = y( k )
373 y( k ) = temp - c( k-1 )*y( k )
374 END IF
375 90 CONTINUE
376 END IF
377*
378* End of SLAGTS
379*
380 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
slagts
subroutine slagts(job, n, a, b, c, d, in, y, tol, info)
SLAGTS solves the system of equations (T-λI)x = y or (T-λI)^Tx = y, where T is a general tridiagonal ...
Definition slagts.f:161