SUBROUTINE DLAGTS( JOB, N, A, B, C, D, IN, Y, TOL, INFO ) * * -- LAPACK auxiliary routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER INFO, JOB, N DOUBLE PRECISION TOL * .. * .. Array Arguments .. INTEGER IN( * ) DOUBLE PRECISION A( * ), B( * ), C( * ), D( * ), Y( * ) * .. * * Purpose * ======= * * DLAGTS may be used to solve one of the systems of equations * * (T - lambda*I)*x = y or (T - lambda*I)'*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. * * Arguments * ========= * * JOB (input) 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)'x = y are to be solved, * but diagonal elements of U are not to be perturbed. * = -2: 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. * * N (input) INTEGER * The order of the matrix T. * * A (input) DOUBLE PRECISION array, dimension (N) * On entry, A must contain the diagonal elements of U as * returned from DLAGTF. * * B (input) DOUBLE PRECISION array, dimension (N-1) * On entry, B must contain the first super-diagonal elements of * U as returned from DLAGTF. * * C (input) DOUBLE PRECISION array, dimension (N-1) * On entry, C must contain the sub-diagonal elements of L as * returned from DLAGTF. * * D (input) DOUBLE PRECISION array, dimension (N-2) * On entry, D must contain the second super-diagonal elements * of U as returned from DLAGTF. * * IN (input) INTEGER array, dimension (N) * On entry, IN must contain details of the matrix P as returned * from DLAGTF. * * Y (input/output) DOUBLE PRECISION array, dimension (N) * On entry, the right hand side vector y. * On exit, Y is overwritten by the solution vector x. * * TOL (input/output) 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. * * INFO (output) 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. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Local Scalars .. INTEGER K DOUBLE PRECISION ABSAK, AK, BIGNUM, EPS, PERT, SFMIN, TEMP * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, SIGN * .. * .. External Functions .. DOUBLE PRECISION DLAMCH EXTERNAL DLAMCH * .. * .. External Subroutines .. EXTERNAL XERBLA * .. * .. Executable Statements .. * INFO = 0 IF( ( ABS( JOB ).GT.2 ) .OR. ( JOB.EQ.0 ) ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DLAGTS', -INFO ) RETURN END IF * IF( N.EQ.0 ) $ RETURN * EPS = DLAMCH( 'Epsilon' ) SFMIN = DLAMCH( 'Safe minimum' ) BIGNUM = ONE / SFMIN * IF( JOB.LT.0 ) THEN IF( TOL.LE.ZERO ) THEN TOL = ABS( A( 1 ) ) IF( N.GT.1 ) $ TOL = MAX( TOL, ABS( A( 2 ) ), ABS( B( 1 ) ) ) DO 10 K = 3, N TOL = MAX( TOL, ABS( A( K ) ), ABS( B( K-1 ) ), $ ABS( D( K-2 ) ) ) 10 CONTINUE TOL = TOL*EPS IF( TOL.EQ.ZERO ) $ TOL = EPS END IF END IF * IF( ABS( JOB ).EQ.1 ) THEN DO 20 K = 2, N IF( IN( K-1 ).EQ.0 ) THEN Y( K ) = Y( K ) - C( K-1 )*Y( K-1 ) ELSE TEMP = Y( K-1 ) Y( K-1 ) = Y( K ) Y( K ) = TEMP - C( K-1 )*Y( K ) END IF 20 CONTINUE IF( JOB.EQ.1 ) THEN DO 30 K = N, 1, -1 IF( K.LE.N-2 ) THEN TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 ) ELSE IF( K.EQ.N-1 ) THEN TEMP = Y( K ) - B( K )*Y( K+1 ) ELSE TEMP = Y( K ) END IF AK = A( K ) ABSAK = ABS( AK ) IF( ABSAK.LT.ONE ) THEN IF( ABSAK.LT.SFMIN ) THEN IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK ) $ THEN INFO = K RETURN ELSE TEMP = TEMP*BIGNUM AK = AK*BIGNUM END IF ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN INFO = K RETURN END IF END IF Y( K ) = TEMP / AK 30 CONTINUE ELSE DO 50 K = N, 1, -1 IF( K.LE.N-2 ) THEN TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 ) ELSE IF( K.EQ.N-1 ) THEN TEMP = Y( K ) - B( K )*Y( K+1 ) ELSE TEMP = Y( K ) END IF AK = A( K ) PERT = SIGN( TOL, AK ) 40 CONTINUE ABSAK = ABS( AK ) IF( ABSAK.LT.ONE ) THEN IF( ABSAK.LT.SFMIN ) THEN IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK ) $ THEN AK = AK + PERT PERT = 2*PERT GO TO 40 ELSE TEMP = TEMP*BIGNUM AK = AK*BIGNUM END IF ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN AK = AK + PERT PERT = 2*PERT GO TO 40 END IF END IF Y( K ) = TEMP / AK 50 CONTINUE END IF ELSE * * Come to here if JOB = 2 or -2 * IF( JOB.EQ.2 ) THEN DO 60 K = 1, N IF( K.GE.3 ) THEN TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 ) ELSE IF( K.EQ.2 ) THEN TEMP = Y( K ) - B( K-1 )*Y( K-1 ) ELSE TEMP = Y( K ) END IF AK = A( K ) ABSAK = ABS( AK ) IF( ABSAK.LT.ONE ) THEN IF( ABSAK.LT.SFMIN ) THEN IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK ) $ THEN INFO = K RETURN ELSE TEMP = TEMP*BIGNUM AK = AK*BIGNUM END IF ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN INFO = K RETURN END IF END IF Y( K ) = TEMP / AK 60 CONTINUE ELSE DO 80 K = 1, N IF( K.GE.3 ) THEN TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 ) ELSE IF( K.EQ.2 ) THEN TEMP = Y( K ) - B( K-1 )*Y( K-1 ) ELSE TEMP = Y( K ) END IF AK = A( K ) PERT = SIGN( TOL, AK ) 70 CONTINUE ABSAK = ABS( AK ) IF( ABSAK.LT.ONE ) THEN IF( ABSAK.LT.SFMIN ) THEN IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK ) $ THEN AK = AK + PERT PERT = 2*PERT GO TO 70 ELSE TEMP = TEMP*BIGNUM AK = AK*BIGNUM END IF ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN AK = AK + PERT PERT = 2*PERT GO TO 70 END IF END IF Y( K ) = TEMP / AK 80 CONTINUE END IF * DO 90 K = N, 2, -1 IF( IN( K-1 ).EQ.0 ) THEN Y( K-1 ) = Y( K-1 ) - C( K-1 )*Y( K ) ELSE TEMP = Y( K-1 ) Y( K-1 ) = Y( K ) Y( K ) = TEMP - C( K-1 )*Y( K ) END IF 90 CONTINUE END IF * * End of DLAGTS * END