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