LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
dlagts.f
Go to the documentation of this file.
1 *> \brief \b DLAGTS solves the system of equations (T-λI)x = y or (T-λI)Tx = y,where T is a general tridiagonal matrix and λ a scalar, using the LU factorization computed by slagtf.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DLAGTS + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlagts.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlagts.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlagts.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLAGTS( JOB, N, A, B, C, D, IN, Y, TOL, INFO )
22 *
23 * .. Scalar Arguments ..
24 * INTEGER INFO, JOB, N
25 * DOUBLE PRECISION TOL
26 * ..
27 * .. Array Arguments ..
28 * INTEGER IN( * )
29 * DOUBLE PRECISION A( * ), B( * ), C( * ), D( * ), Y( * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> DLAGTS 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 DLAGTF. 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 DLAGTS 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 DOUBLE PRECISION array, dimension (N)
83 *> On entry, A must contain the diagonal elements of U as
84 *> returned from DLAGTF.
85 *> \endverbatim
86 *>
87 *> \param[in] B
88 *> \verbatim
89 *> B is DOUBLE PRECISION array, dimension (N-1)
90 *> On entry, B must contain the first super-diagonal elements of
91 *> U as returned from DLAGTF.
92 *> \endverbatim
93 *>
94 *> \param[in] C
95 *> \verbatim
96 *> C is DOUBLE PRECISION array, dimension (N-1)
97 *> On entry, C must contain the sub-diagonal elements of L as
98 *> returned from DLAGTF.
99 *> \endverbatim
100 *>
101 *> \param[in] D
102 *> \verbatim
103 *> D is DOUBLE PRECISION array, dimension (N-2)
104 *> On entry, D must contain the second super-diagonal elements
105 *> of U as returned from DLAGTF.
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 DLAGTF.
113 *> \endverbatim
114 *>
115 *> \param[in,out] Y
116 *> \verbatim
117 *> Y is DOUBLE PRECISION 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 DOUBLE PRECISION
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 OTHERauxiliary
158 *
159 * =====================================================================
160  SUBROUTINE dlagts( 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  DOUBLE PRECISION TOL
169 * ..
170 * .. Array Arguments ..
171  INTEGER IN( * )
172  DOUBLE PRECISION A( * ), B( * ), C( * ), D( * ), Y( * )
173 * ..
174 *
175 * =====================================================================
176 *
177 * .. Parameters ..
178  DOUBLE PRECISION ONE, ZERO
179  parameter( one = 1.0d+0, zero = 0.0d+0 )
180 * ..
181 * .. Local Scalars ..
182  INTEGER K
183  DOUBLE PRECISION ABSAK, AK, BIGNUM, EPS, PERT, SFMIN, TEMP
184 * ..
185 * .. Intrinsic Functions ..
186  INTRINSIC abs, max, sign
187 * ..
188 * .. External Functions ..
189  DOUBLE PRECISION DLAMCH
190  EXTERNAL dlamch
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( 'DLAGTS', -info )
205  RETURN
206  END IF
207 *
208  IF( n.EQ.0 )
209  $ RETURN
210 *
211  eps = dlamch( 'Epsilon' )
212  sfmin = dlamch( '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 DLAGTS
379 *
380  END
subroutine dlagts(JOB, N, A, B, C, D, IN, Y, TOL, INFO)
DLAGTS solves the system of equations (T-λI)x = y or (T-λI)Tx = y,where T is a general tridiagonal ma...
Definition: dlagts.f:161
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60