LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ dlagts()

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 < 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 > 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
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 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.

Definition at line 160 of file dlagts.f.

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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
Here is the call graph for this function:
Here is the caller graph for this function: