diff options
Diffstat (limited to 'SRC/slagts.f')
| -rw-r--r-- | SRC/slagts.f | 304 |
1 files changed, 304 insertions, 0 deletions
diff --git a/SRC/slagts.f b/SRC/slagts.f new file mode 100644 index 00000000..e2f43bee --- /dev/null +++ b/SRC/slagts.f @@ -0,0 +1,304 @@ + SUBROUTINE SLAGTS( 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 + REAL TOL +* .. +* .. Array Arguments .. + INTEGER IN( * ) + REAL A( * ), B( * ), C( * ), D( * ), Y( * ) +* .. +* +* Purpose +* ======= +* +* SLAGTS 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 SLAGTF. 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 SLAGTS 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) REAL array, dimension (N) +* On entry, A must contain the diagonal elements of U as +* returned from SLAGTF. +* +* B (input) REAL array, dimension (N-1) +* On entry, B must contain the first super-diagonal elements of +* U as returned from SLAGTF. +* +* C (input) REAL array, dimension (N-1) +* On entry, C must contain the sub-diagonal elements of L as +* returned from SLAGTF. +* +* D (input) REAL array, dimension (N-2) +* On entry, D must contain the second super-diagonal elements +* of U as returned from SLAGTF. +* +* IN (input) INTEGER array, dimension (N) +* On entry, IN must contain details of the matrix P as returned +* from SLAGTF. +* +* Y (input/output) REAL array, dimension (N) +* On entry, the right hand side vector y. +* On exit, Y is overwritten by the solution vector x. +* +* TOL (input/output) REAL +* 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 .. + REAL ONE, ZERO + PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) +* .. +* .. Local Scalars .. + INTEGER K + REAL ABSAK, AK, BIGNUM, EPS, PERT, SFMIN, TEMP +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, MAX, SIGN +* .. +* .. External Functions .. + REAL SLAMCH + EXTERNAL SLAMCH +* .. +* .. 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( 'SLAGTS', -INFO ) + RETURN + END IF +* + IF( N.EQ.0 ) + $ RETURN +* + EPS = SLAMCH( 'Epsilon' ) + SFMIN = SLAMCH( '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 SLAGTS +* + END |