LAPACK: SRC/dlagts.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DLAGTS( JOB, N, A, B, C, D, IN, Y, TOL, INFO )
00002 *
00003 *  -- LAPACK auxiliary routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            INFO, JOB, N
00010       DOUBLE PRECISION   TOL
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            IN( * )
00014       DOUBLE PRECISION   A( * ), B( * ), C( * ), D( * ), Y( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  DLAGTS may be used to solve one of the systems of equations
00021 *
00022 *     (T - lambda*I)*x = y   or   (T - lambda*I)'*x = y,
00023 *
00024 *  where T is an n by n tridiagonal matrix, for x, following the
00025 *  factorization of (T - lambda*I) as
00026 *
00027 *     (T - lambda*I) = P*L*U ,
00028 *
00029 *  by routine DLAGTF. The choice of equation to be solved is
00030 *  controlled by the argument JOB, and in each case there is an option
00031 *  to perturb zero or very small diagonal elements of U, this option
00032 *  being intended for use in applications such as inverse iteration.
00033 *
00034 *  Arguments
00035 *  =========
00036 *
00037 *  JOB     (input) INTEGER
00038 *          Specifies the job to be performed by DLAGTS as follows:
00039 *          =  1: The equations  (T - lambda*I)x = y  are to be solved,
00040 *                but diagonal elements of U are not to be perturbed.
00041 *          = -1: The equations  (T - lambda*I)x = y  are to be solved
00042 *                and, if overflow would otherwise occur, the diagonal
00043 *                elements of U are to be perturbed. See argument TOL
00044 *                below.
00045 *          =  2: The equations  (T - lambda*I)'x = y  are to be solved,
00046 *                but diagonal elements of U are not to be perturbed.
00047 *          = -2: The equations  (T - lambda*I)'x = y  are to be solved
00048 *                and, if overflow would otherwise occur, the diagonal
00049 *                elements of U are to be perturbed. See argument TOL
00050 *                below.
00051 *
00052 *  N       (input) INTEGER
00053 *          The order of the matrix T.
00054 *
00055 *  A       (input) DOUBLE PRECISION array, dimension (N)
00056 *          On entry, A must contain the diagonal elements of U as
00057 *          returned from DLAGTF.
00058 *
00059 *  B       (input) DOUBLE PRECISION array, dimension (N-1)
00060 *          On entry, B must contain the first super-diagonal elements of
00061 *          U as returned from DLAGTF.
00062 *
00063 *  C       (input) DOUBLE PRECISION array, dimension (N-1)
00064 *          On entry, C must contain the sub-diagonal elements of L as
00065 *          returned from DLAGTF.
00066 *
00067 *  D       (input) DOUBLE PRECISION array, dimension (N-2)
00068 *          On entry, D must contain the second super-diagonal elements
00069 *          of U as returned from DLAGTF.
00070 *
00071 *  IN      (input) INTEGER array, dimension (N)
00072 *          On entry, IN must contain details of the matrix P as returned
00073 *          from DLAGTF.
00074 *
00075 *  Y       (input/output) DOUBLE PRECISION array, dimension (N)
00076 *          On entry, the right hand side vector y.
00077 *          On exit, Y is overwritten by the solution vector x.
00078 *
00079 *  TOL     (input/output) DOUBLE PRECISION
00080 *          On entry, with  JOB .lt. 0, TOL should be the minimum
00081 *          perturbation to be made to very small diagonal elements of U.
00082 *          TOL should normally be chosen as about eps*norm(U), where eps
00083 *          is the relative machine precision, but if TOL is supplied as
00084 *          non-positive, then it is reset to eps*max( abs( u(i,j) ) ).
00085 *          If  JOB .gt. 0  then TOL is not referenced.
00086 *
00087 *          On exit, TOL is changed as described above, only if TOL is
00088 *          non-positive on entry. Otherwise TOL is unchanged.
00089 *
00090 *  INFO    (output) INTEGER
00091 *          = 0   : successful exit
00092 *          .lt. 0: if INFO = -i, the i-th argument had an illegal value
00093 *          .gt. 0: overflow would occur when computing the INFO(th)
00094 *                  element of the solution vector x. This can only occur
00095 *                  when JOB is supplied as positive and either means
00096 *                  that a diagonal element of U is very small, or that
00097 *                  the elements of the right-hand side vector y are very
00098 *                  large.
00099 *
00100 *  =====================================================================
00101 *
00102 *     .. Parameters ..
00103       DOUBLE PRECISION   ONE, ZERO
00104       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00105 *     ..
00106 *     .. Local Scalars ..
00107       INTEGER            K
00108       DOUBLE PRECISION   ABSAK, AK, BIGNUM, EPS, PERT, SFMIN, TEMP
00109 *     ..
00110 *     .. Intrinsic Functions ..
00111       INTRINSIC          ABS, MAX, SIGN
00112 *     ..
00113 *     .. External Functions ..
00114       DOUBLE PRECISION   DLAMCH
00115       EXTERNAL           DLAMCH
00116 *     ..
00117 *     .. External Subroutines ..
00118       EXTERNAL           XERBLA
00119 *     ..
00120 *     .. Executable Statements ..
00121 *
00122       INFO = 0
00123       IF( ( ABS( JOB ).GT.2 ) .OR. ( JOB.EQ.0 ) ) THEN
00124          INFO = -1
00125       ELSE IF( N.LT.0 ) THEN
00126          INFO = -2
00127       END IF
00128       IF( INFO.NE.0 ) THEN
00129          CALL XERBLA( 'DLAGTS', -INFO )
00130          RETURN
00131       END IF
00132 *
00133       IF( N.EQ.0 )
00134      $   RETURN
00135 *
00136       EPS = DLAMCH( 'Epsilon' )
00137       SFMIN = DLAMCH( 'Safe minimum' )
00138       BIGNUM = ONE / SFMIN
00139 *
00140       IF( JOB.LT.0 ) THEN
00141          IF( TOL.LE.ZERO ) THEN
00142             TOL = ABS( A( 1 ) )
00143             IF( N.GT.1 )
00144      $         TOL = MAX( TOL, ABS( A( 2 ) ), ABS( B( 1 ) ) )
00145             DO 10 K = 3, N
00146                TOL = MAX( TOL, ABS( A( K ) ), ABS( B( K-1 ) ),
00147      $               ABS( D( K-2 ) ) )
00148    10       CONTINUE
00149             TOL = TOL*EPS
00150             IF( TOL.EQ.ZERO )
00151      $         TOL = EPS
00152          END IF
00153       END IF
00154 *
00155       IF( ABS( JOB ).EQ.1 ) THEN
00156          DO 20 K = 2, N
00157             IF( IN( K-1 ).EQ.0 ) THEN
00158                Y( K ) = Y( K ) - C( K-1 )*Y( K-1 )
00159             ELSE
00160                TEMP = Y( K-1 )
00161                Y( K-1 ) = Y( K )
00162                Y( K ) = TEMP - C( K-1 )*Y( K )
00163             END IF
00164    20    CONTINUE
00165          IF( JOB.EQ.1 ) THEN
00166             DO 30 K = N, 1, -1
00167                IF( K.LE.N-2 ) THEN
00168                   TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 )
00169                ELSE IF( K.EQ.N-1 ) THEN
00170                   TEMP = Y( K ) - B( K )*Y( K+1 )
00171                ELSE
00172                   TEMP = Y( K )
00173                END IF
00174                AK = A( K )
00175                ABSAK = ABS( AK )
00176                IF( ABSAK.LT.ONE ) THEN
00177                   IF( ABSAK.LT.SFMIN ) THEN
00178                      IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
00179      $                    THEN
00180                         INFO = K
00181                         RETURN
00182                      ELSE
00183                         TEMP = TEMP*BIGNUM
00184                         AK = AK*BIGNUM
00185                      END IF
00186                   ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
00187                      INFO = K
00188                      RETURN
00189                   END IF
00190                END IF
00191                Y( K ) = TEMP / AK
00192    30       CONTINUE
00193          ELSE
00194             DO 50 K = N, 1, -1
00195                IF( K.LE.N-2 ) THEN
00196                   TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 )
00197                ELSE IF( K.EQ.N-1 ) THEN
00198                   TEMP = Y( K ) - B( K )*Y( K+1 )
00199                ELSE
00200                   TEMP = Y( K )
00201                END IF
00202                AK = A( K )
00203                PERT = SIGN( TOL, AK )
00204    40          CONTINUE
00205                ABSAK = ABS( AK )
00206                IF( ABSAK.LT.ONE ) THEN
00207                   IF( ABSAK.LT.SFMIN ) THEN
00208                      IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
00209      $                    THEN
00210                         AK = AK + PERT
00211                         PERT = 2*PERT
00212                         GO TO 40
00213                      ELSE
00214                         TEMP = TEMP*BIGNUM
00215                         AK = AK*BIGNUM
00216                      END IF
00217                   ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
00218                      AK = AK + PERT
00219                      PERT = 2*PERT
00220                      GO TO 40
00221                   END IF
00222                END IF
00223                Y( K ) = TEMP / AK
00224    50       CONTINUE
00225          END IF
00226       ELSE
00227 *
00228 *        Come to here if  JOB = 2 or -2
00229 *
00230          IF( JOB.EQ.2 ) THEN
00231             DO 60 K = 1, N
00232                IF( K.GE.3 ) THEN
00233                   TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 )
00234                ELSE IF( K.EQ.2 ) THEN
00235                   TEMP = Y( K ) - B( K-1 )*Y( K-1 )
00236                ELSE
00237                   TEMP = Y( K )
00238                END IF
00239                AK = A( K )
00240                ABSAK = ABS( AK )
00241                IF( ABSAK.LT.ONE ) THEN
00242                   IF( ABSAK.LT.SFMIN ) THEN
00243                      IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
00244      $                    THEN
00245                         INFO = K
00246                         RETURN
00247                      ELSE
00248                         TEMP = TEMP*BIGNUM
00249                         AK = AK*BIGNUM
00250                      END IF
00251                   ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
00252                      INFO = K
00253                      RETURN
00254                   END IF
00255                END IF
00256                Y( K ) = TEMP / AK
00257    60       CONTINUE
00258          ELSE
00259             DO 80 K = 1, N
00260                IF( K.GE.3 ) THEN
00261                   TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 )
00262                ELSE IF( K.EQ.2 ) THEN
00263                   TEMP = Y( K ) - B( K-1 )*Y( K-1 )
00264                ELSE
00265                   TEMP = Y( K )
00266                END IF
00267                AK = A( K )
00268                PERT = SIGN( TOL, AK )
00269    70          CONTINUE
00270                ABSAK = ABS( AK )
00271                IF( ABSAK.LT.ONE ) THEN
00272                   IF( ABSAK.LT.SFMIN ) THEN
00273                      IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
00274      $                    THEN
00275                         AK = AK + PERT
00276                         PERT = 2*PERT
00277                         GO TO 70
00278                      ELSE
00279                         TEMP = TEMP*BIGNUM
00280                         AK = AK*BIGNUM
00281                      END IF
00282                   ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
00283                      AK = AK + PERT
00284                      PERT = 2*PERT
00285                      GO TO 70
00286                   END IF
00287                END IF
00288                Y( K ) = TEMP / AK
00289    80       CONTINUE
00290          END IF
00291 *
00292          DO 90 K = N, 2, -1
00293             IF( IN( K-1 ).EQ.0 ) THEN
00294                Y( K-1 ) = Y( K-1 ) - C( K-1 )*Y( K )
00295             ELSE
00296                TEMP = Y( K-1 )
00297                Y( K-1 ) = Y( K )
00298                Y( K ) = TEMP - C( K-1 )*Y( K )
00299             END IF
00300    90    CONTINUE
00301       END IF
00302 *
00303 *     End of DLAGTS
00304 *
00305       END