LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 .lt. 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 .gt. 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 *> .lt. 0: if INFO = -i, the i-th argument had an illegal value
141 *> .gt. 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 *> \date September 2012
158 *
159 *> \ingroup auxOTHERauxiliary
160 *
161 * =====================================================================
162  SUBROUTINE dlagts( JOB, N, A, B, C, D, IN, Y, TOL, INFO )
163 *
164 * -- LAPACK auxiliary routine (version 3.4.2) --
165 * -- LAPACK is a software package provided by Univ. of Tennessee, --
166 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
167 * September 2012
168 *
169 * .. Scalar Arguments ..
170  INTEGER info, job, n
171  DOUBLE PRECISION tol
172 * ..
173 * .. Array Arguments ..
174  INTEGER in( * )
175  DOUBLE PRECISION a( * ), b( * ), c( * ), d( * ), y( * )
176 * ..
177 *
178 * =====================================================================
179 *
180 * .. Parameters ..
181  DOUBLE PRECISION one, zero
182  parameter( one = 1.0d+0, zero = 0.0d+0 )
183 * ..
184 * .. Local Scalars ..
185  INTEGER k
186  DOUBLE PRECISION absak, ak, bignum, eps, pert, sfmin, temp
187 * ..
188 * .. Intrinsic Functions ..
189  INTRINSIC abs, max, sign
190 * ..
191 * .. External Functions ..
192  DOUBLE PRECISION dlamch
193  EXTERNAL dlamch
194 * ..
195 * .. External Subroutines ..
196  EXTERNAL xerbla
197 * ..
198 * .. Executable Statements ..
199 *
200  info = 0
201  IF( ( abs( job ).GT.2 ) .OR. ( job.EQ.0 ) ) THEN
202  info = -1
203  ELSE IF( n.LT.0 ) THEN
204  info = -2
205  END IF
206  IF( info.NE.0 ) THEN
207  CALL xerbla( 'DLAGTS', -info )
208  return
209  END IF
210 *
211  IF( n.EQ.0 )
212  $ return
213 *
214  eps = dlamch( 'Epsilon' )
215  sfmin = dlamch( 'Safe minimum' )
216  bignum = one / sfmin
217 *
218  IF( job.LT.0 ) THEN
219  IF( tol.LE.zero ) THEN
220  tol = abs( a( 1 ) )
221  IF( n.GT.1 )
222  $ tol = max( tol, abs( a( 2 ) ), abs( b( 1 ) ) )
223  DO 10 k = 3, n
224  tol = max( tol, abs( a( k ) ), abs( b( k-1 ) ),
225  $ abs( d( k-2 ) ) )
226  10 continue
227  tol = tol*eps
228  IF( tol.EQ.zero )
229  $ tol = eps
230  END IF
231  END IF
232 *
233  IF( abs( job ).EQ.1 ) THEN
234  DO 20 k = 2, n
235  IF( in( k-1 ).EQ.0 ) THEN
236  y( k ) = y( k ) - c( k-1 )*y( k-1 )
237  ELSE
238  temp = y( k-1 )
239  y( k-1 ) = y( k )
240  y( k ) = temp - c( k-1 )*y( k )
241  END IF
242  20 continue
243  IF( job.EQ.1 ) THEN
244  DO 30 k = n, 1, -1
245  IF( k.LE.n-2 ) THEN
246  temp = y( k ) - b( k )*y( k+1 ) - d( k )*y( k+2 )
247  ELSE IF( k.EQ.n-1 ) THEN
248  temp = y( k ) - b( k )*y( k+1 )
249  ELSE
250  temp = y( k )
251  END IF
252  ak = a( k )
253  absak = abs( ak )
254  IF( absak.LT.one ) THEN
255  IF( absak.LT.sfmin ) THEN
256  IF( absak.EQ.zero .OR. abs( temp )*sfmin.GT.absak )
257  $ THEN
258  info = k
259  return
260  ELSE
261  temp = temp*bignum
262  ak = ak*bignum
263  END IF
264  ELSE IF( abs( temp ).GT.absak*bignum ) THEN
265  info = k
266  return
267  END IF
268  END IF
269  y( k ) = temp / ak
270  30 continue
271  ELSE
272  DO 50 k = n, 1, -1
273  IF( k.LE.n-2 ) THEN
274  temp = y( k ) - b( k )*y( k+1 ) - d( k )*y( k+2 )
275  ELSE IF( k.EQ.n-1 ) THEN
276  temp = y( k ) - b( k )*y( k+1 )
277  ELSE
278  temp = y( k )
279  END IF
280  ak = a( k )
281  pert = sign( tol, ak )
282  40 continue
283  absak = abs( ak )
284  IF( absak.LT.one ) THEN
285  IF( absak.LT.sfmin ) THEN
286  IF( absak.EQ.zero .OR. abs( temp )*sfmin.GT.absak )
287  $ THEN
288  ak = ak + pert
289  pert = 2*pert
290  go to 40
291  ELSE
292  temp = temp*bignum
293  ak = ak*bignum
294  END IF
295  ELSE IF( abs( temp ).GT.absak*bignum ) THEN
296  ak = ak + pert
297  pert = 2*pert
298  go to 40
299  END IF
300  END IF
301  y( k ) = temp / ak
302  50 continue
303  END IF
304  ELSE
305 *
306 * Come to here if JOB = 2 or -2
307 *
308  IF( job.EQ.2 ) THEN
309  DO 60 k = 1, n
310  IF( k.GE.3 ) THEN
311  temp = y( k ) - b( k-1 )*y( k-1 ) - d( k-2 )*y( k-2 )
312  ELSE IF( k.EQ.2 ) THEN
313  temp = y( k ) - b( k-1 )*y( k-1 )
314  ELSE
315  temp = y( k )
316  END IF
317  ak = a( k )
318  absak = abs( ak )
319  IF( absak.LT.one ) THEN
320  IF( absak.LT.sfmin ) THEN
321  IF( absak.EQ.zero .OR. abs( temp )*sfmin.GT.absak )
322  $ THEN
323  info = k
324  return
325  ELSE
326  temp = temp*bignum
327  ak = ak*bignum
328  END IF
329  ELSE IF( abs( temp ).GT.absak*bignum ) THEN
330  info = k
331  return
332  END IF
333  END IF
334  y( k ) = temp / ak
335  60 continue
336  ELSE
337  DO 80 k = 1, n
338  IF( k.GE.3 ) THEN
339  temp = y( k ) - b( k-1 )*y( k-1 ) - d( k-2 )*y( k-2 )
340  ELSE IF( k.EQ.2 ) THEN
341  temp = y( k ) - b( k-1 )*y( k-1 )
342  ELSE
343  temp = y( k )
344  END IF
345  ak = a( k )
346  pert = sign( tol, ak )
347  70 continue
348  absak = abs( ak )
349  IF( absak.LT.one ) THEN
350  IF( absak.LT.sfmin ) THEN
351  IF( absak.EQ.zero .OR. abs( temp )*sfmin.GT.absak )
352  $ THEN
353  ak = ak + pert
354  pert = 2*pert
355  go to 70
356  ELSE
357  temp = temp*bignum
358  ak = ak*bignum
359  END IF
360  ELSE IF( abs( temp ).GT.absak*bignum ) THEN
361  ak = ak + pert
362  pert = 2*pert
363  go to 70
364  END IF
365  END IF
366  y( k ) = temp / ak
367  80 continue
368  END IF
369 *
370  DO 90 k = n, 2, -1
371  IF( in( k-1 ).EQ.0 ) THEN
372  y( k-1 ) = y( k-1 ) - c( k-1 )*y( k )
373  ELSE
374  temp = y( k-1 )
375  y( k-1 ) = y( k )
376  y( k ) = temp - c( k-1 )*y( k )
377  END IF
378  90 continue
379  END IF
380 *
381 * End of DLAGTS
382 *
383  END