LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dlagts ( integer  JOB,
integer  N,
double precision, dimension( * )  A,
double precision, dimension( * )  B,
double precision, dimension( * )  C,
double precision, dimension( * )  D,
integer, dimension( * )  IN,
double precision, dimension( * )  Y,
double precision  TOL,
integer  INFO 
)

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.

Download DLAGTS + dependencies [TGZ] [ZIP] [TXT]

Purpose: 
 DLAGTS may be used to solve one of the systems of equations

    (T - lambda*I)*x = y   or   (T - lambda*I)**T*x = y,

 where T is an n by n tridiagonal matrix, for x, following the
 factorization of (T - lambda*I) as

    (T - lambda*I) = P*L*U ,

 by routine DLAGTF. The choice of equation to be solved is
 controlled by the argument JOB, and in each case there is an option
 to perturb zero or very small diagonal elements of U, this option
 being intended for use in applications such as inverse iteration.
Parameters
[in]JOB
          JOB is INTEGER
          Specifies the job to be performed by DLAGTS as follows:
          =  1: The equations  (T - lambda*I)x = y  are to be solved,
                but diagonal elements of U are not to be perturbed.
          = -1: The equations  (T - lambda*I)x = y  are to be solved
                and, if overflow would otherwise occur, the diagonal
                elements of U are to be perturbed. See argument TOL
                below.
          =  2: The equations  (T - lambda*I)**Tx = y  are to be solved,
                but diagonal elements of U are not to be perturbed.
          = -2: The equations  (T - lambda*I)**Tx = y  are to be solved
                and, if overflow would otherwise occur, the diagonal
                elements of U are to be perturbed. See argument TOL
                below.
[in]N
          N is INTEGER
          The order of the matrix T.
[in]A
          A is DOUBLE PRECISION array, dimension (N)
          On entry, A must contain the diagonal elements of U as
          returned from DLAGTF.
[in]B
          B is DOUBLE PRECISION array, dimension (N-1)
          On entry, B must contain the first super-diagonal elements of
          U as returned from DLAGTF.
[in]C
          C is DOUBLE PRECISION array, dimension (N-1)
          On entry, C must contain the sub-diagonal elements of L as
          returned from DLAGTF.
[in]D
          D is DOUBLE PRECISION array, dimension (N-2)
          On entry, D must contain the second super-diagonal elements
          of U as returned from DLAGTF.
[in]IN
          IN is INTEGER array, dimension (N)
          On entry, IN must contain details of the matrix P as returned
          from DLAGTF.
[in,out]Y
          Y is DOUBLE PRECISION array, dimension (N)
          On entry, the right hand side vector y.
          On exit, Y is overwritten by the solution vector x.
[in,out]TOL
          TOL is DOUBLE PRECISION
          On entry, with  JOB .lt. 0, TOL should be the minimum
          perturbation to be made to very small diagonal elements of U.
          TOL should normally be chosen as about eps*norm(U), where eps
          is the relative machine precision, but if TOL is supplied as
          non-positive, then it is reset to eps*max( abs( u(i,j) ) ).
          If  JOB .gt. 0  then TOL is not referenced.

          On exit, TOL is changed as described above, only if TOL is
          non-positive on entry. Otherwise TOL is unchanged.
[out]INFO
          INFO is INTEGER
          = 0   : successful exit
          .lt. 0: if INFO = -i, the i-th argument had an illegal value
          .gt. 0: overflow would occur when computing the INFO(th)
                  element of the solution vector x. This can only occur
                  when JOB is supplied as positive and either means
                  that a diagonal element of U is very small, or that
                  the elements of the right-hand side vector y are very
                  large.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 163 of file dlagts.f.

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 *
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62

Here is the call graph for this function:

Here is the caller graph for this function: