src/LAPACK/dlagts.f

   1       SUBROUTINE DLAGTS( JOB, N, A, B, C, D, IN, Y, TOL, INFO )
   2 *
   3 *  -- LAPACK auxiliary routine (version 3.1) --
   4 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
   5 *     November 2006
   6 *
   7 *     .. Scalar Arguments ..
   8       INTEGER            INFO, JOB, N
   9       DOUBLE PRECISION   TOL
  10 *     ..
  11 *     .. Array Arguments ..
  12       INTEGER            IN( * )
  13       DOUBLE PRECISION   A( * ), B( * ), C( * ), D( * ), Y( * )
  14 *     ..
  15 *
  16 *  Purpose
  17 *  =======
  18 *
  19 *  DLAGTS may be used to solve one of the systems of equations
  20 *
  21 *     (T - lambda*I)*x = y   or   (T - lambda*I)'*x = y,
  22 *
  23 *  where T is an n by n tridiagonal matrix, for x, following the
  24 *  factorization of (T - lambda*I) as
  25 *
  26 *     (T - lambda*I) = P*L*U ,
  27 *
  28 *  by routine DLAGTF. The choice of equation to be solved is
  29 *  controlled by the argument JOB, and in each case there is an option
  30 *  to perturb zero or very small diagonal elements of U, this option
  31 *  being intended for use in applications such as inverse iteration.
  32 *
  33 *  Arguments
  34 *  =========
  35 *
  36 *  JOB     (input) INTEGER
  37 *          Specifies the job to be performed by DLAGTS as follows:
  38 *          =  1: The equations  (T - lambda*I)x = y  are to be solved,
  39 *                but diagonal elements of U are not to be perturbed.
  40 *          = -1: The equations  (T - lambda*I)x = y  are to be solved
  41 *                and, if overflow would otherwise occur, the diagonal
  42 *                elements of U are to be perturbed. See argument TOL
  43 *                below.
  44 *          =  2: The equations  (T - lambda*I)'x = y  are to be solved,
  45 *                but diagonal elements of U are not to be perturbed.
  46 *          = -2: The equations  (T - lambda*I)'x = y  are to be solved
  47 *                and, if overflow would otherwise occur, the diagonal
  48 *                elements of U are to be perturbed. See argument TOL
  49 *                below.
  50 *
  51 *  N       (input) INTEGER
  52 *          The order of the matrix T.
  53 *
  54 *  A       (input) DOUBLE PRECISION array, dimension (N)
  55 *          On entry, A must contain the diagonal elements of U as
  56 *          returned from DLAGTF.
  57 *
  58 *  B       (input) DOUBLE PRECISION array, dimension (N-1)
  59 *          On entry, B must contain the first super-diagonal elements of
  60 *          U as returned from DLAGTF.
  61 *
  62 *  C       (input) DOUBLE PRECISION array, dimension (N-1)
  63 *          On entry, C must contain the sub-diagonal elements of L as
  64 *          returned from DLAGTF.
  65 *
  66 *  D       (input) DOUBLE PRECISION array, dimension (N-2)
  67 *          On entry, D must contain the second super-diagonal elements
  68 *          of U as returned from DLAGTF.
  69 *
  70 *  IN      (input) INTEGER array, dimension (N)
  71 *          On entry, IN must contain details of the matrix P as returned
  72 *          from DLAGTF.
  73 *
  74 *  Y       (input/output) DOUBLE PRECISION array, dimension (N)
  75 *          On entry, the right hand side vector y.
  76 *          On exit, Y is overwritten by the solution vector x.
  77 *
  78 *  TOL     (input/output) DOUBLE PRECISION
  79 *          On entry, with  JOB .lt. 0, TOL should be the minimum
  80 *          perturbation to be made to very small diagonal elements of U.
  81 *          TOL should normally be chosen as about eps*norm(U), where eps
  82 *          is the relative machine precision, but if TOL is supplied as
  83 *          non-positive, then it is reset to eps*max( abs( u(i,j) ) ).
  84 *          If  JOB .gt. 0  then TOL is not referenced.
  85 *
  86 *          On exit, TOL is changed as described above, only if TOL is
  87 *          non-positive on entry. Otherwise TOL is unchanged.
  88 *
  89 *  INFO    (output) INTEGER
  90 *          = 0   : successful exit
  91 *          .lt. 0: if INFO = -i, the i-th argument had an illegal value
  92 *          .gt. 0: overflow would occur when computing the INFO(th)
  93 *                  element of the solution vector x. This can only occur
  94 *                  when JOB is supplied as positive and either means
  95 *                  that a diagonal element of U is very small, or that
  96 *                  the elements of the right-hand side vector y are very
  97 *                  large.
  98 *
  99 *  =====================================================================
 100 *
 101 *     .. Parameters ..
 102       DOUBLE PRECISION   ONE, ZERO
 103       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 104 *     ..
 105 *     .. Local Scalars ..
 106       INTEGER            K
 107       DOUBLE PRECISION   ABSAK, AK, BIGNUM, EPS, PERT, SFMIN, TEMP
 108 *     ..
 109 *     .. Intrinsic Functions ..
 110       INTRINSIC          ABS, MAX, SIGN
 111 *     ..
 112 *     .. External Functions ..
 113       DOUBLE PRECISION   DLAMCH
 114       EXTERNAL           DLAMCH
 115 *     ..
 116 *     .. External Subroutines ..
 117       EXTERNAL           XERBLA
 118 *     ..
 119 *     .. Executable Statements ..
 120 *
 121       INFO = 0
 122       IF( ( ABS( JOB ).GT.2 ) .OR. ( JOB.EQ.0 ) ) THEN
 123          INFO = -1
 124       ELSE IF( N.LT.0 ) THEN
 125          INFO = -2
 126       END IF
 127       IF( INFO.NE.0 ) THEN
 128          CALL XERBLA( 'DLAGTS', -INFO )
 129          RETURN
 130       END IF
 131 *
 132       IF( N.EQ.0 )
 133      $   RETURN
 134 *
 135       EPS = DLAMCH( 'Epsilon' )
 136       SFMIN = DLAMCH( 'Safe minimum' )
 137       BIGNUM = ONE / SFMIN
 138 *
 139       IF( JOB.LT.0 ) THEN
 140          IF( TOL.LE.ZERO ) THEN
 141             TOL = ABS( A( 1 ) )
 142             IF( N.GT.1 )
 143      $         TOL = MAX( TOL, ABS( A( 2 ) ), ABS( B( 1 ) ) )
 144             DO 10 K = 3, N
 145                TOL = MAX( TOL, ABS( A( K ) ), ABS( B( K-1 ) ),
 146      $               ABS( D( K-2 ) ) )
 147    10       CONTINUE
 148             TOL = TOL*EPS
 149             IF( TOL.EQ.ZERO )
 150      $         TOL = EPS
 151          END IF
 152       END IF
 153 *
 154       IF( ABS( JOB ).EQ.1 ) THEN
 155          DO 20 K = 2, N
 156             IF( IN( K-1 ).EQ.0 ) THEN
 157                Y( K ) = Y( K ) - C( K-1 )*Y( K-1 )
 158             ELSE
 159                TEMP = Y( K-1 )
 160                Y( K-1 ) = Y( K )
 161                Y( K ) = TEMP - C( K-1 )*Y( K )
 162             END IF
 163    20    CONTINUE
 164          IF( JOB.EQ.1 ) THEN
 165             DO 30 K = N, 1, -1
 166                IF( K.LE.N-2 ) THEN
 167                   TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 )
 168                ELSE IF( K.EQ.N-1 ) THEN
 169                   TEMP = Y( K ) - B( K )*Y( K+1 )
 170                ELSE
 171                   TEMP = Y( K )
 172                END IF
 173                AK = A( K )
 174                ABSAK = ABS( AK )
 175                IF( ABSAK.LT.ONE ) THEN
 176                   IF( ABSAK.LT.SFMIN ) THEN
 177                      IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
 178      $                    THEN
 179                         INFO = K
 180                         RETURN
 181                      ELSE
 182                         TEMP = TEMP*BIGNUM
 183                         AK = AK*BIGNUM
 184                      END IF
 185                   ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
 186                      INFO = K
 187                      RETURN
 188                   END IF
 189                END IF
 190                Y( K ) = TEMP / AK
 191    30       CONTINUE
 192          ELSE
 193             DO 50 K = N, 1, -1
 194                IF( K.LE.N-2 ) THEN
 195                   TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 )
 196                ELSE IF( K.EQ.N-1 ) THEN
 197                   TEMP = Y( K ) - B( K )*Y( K+1 )
 198                ELSE
 199                   TEMP = Y( K )
 200                END IF
 201                AK = A( K )
 202                PERT = SIGN( TOL, AK )
 203    40          CONTINUE
 204                ABSAK = ABS( AK )
 205                IF( ABSAK.LT.ONE ) THEN
 206                   IF( ABSAK.LT.SFMIN ) THEN
 207                      IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
 208      $                    THEN
 209                         AK = AK + PERT
 210                         PERT = 2*PERT
 211                         GO TO 40
 212                      ELSE
 213                         TEMP = TEMP*BIGNUM
 214                         AK = AK*BIGNUM
 215                      END IF
 216                   ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
 217                      AK = AK + PERT
 218                      PERT = 2*PERT
 219                      GO TO 40
 220                   END IF
 221                END IF
 222                Y( K ) = TEMP / AK
 223    50       CONTINUE
 224          END IF
 225       ELSE
 226 *
 227 *        Come to here if  JOB = 2 or -2
 228 *
 229          IF( JOB.EQ.2 ) THEN
 230             DO 60 K = 1, N
 231                IF( K.GE.3 ) THEN
 232                   TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 )
 233                ELSE IF( K.EQ.2 ) THEN
 234                   TEMP = Y( K ) - B( K-1 )*Y( K-1 )
 235                ELSE
 236                   TEMP = Y( K )
 237                END IF
 238                AK = A( K )
 239                ABSAK = ABS( AK )
 240                IF( ABSAK.LT.ONE ) THEN
 241                   IF( ABSAK.LT.SFMIN ) THEN
 242                      IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
 243      $                    THEN
 244                         INFO = K
 245                         RETURN
 246                      ELSE
 247                         TEMP = TEMP*BIGNUM
 248                         AK = AK*BIGNUM
 249                      END IF
 250                   ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
 251                      INFO = K
 252                      RETURN
 253                   END IF
 254                END IF
 255                Y( K ) = TEMP / AK
 256    60       CONTINUE
 257          ELSE
 258             DO 80 K = 1, N
 259                IF( K.GE.3 ) THEN
 260                   TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 )
 261                ELSE IF( K.EQ.2 ) THEN
 262                   TEMP = Y( K ) - B( K-1 )*Y( K-1 )
 263                ELSE
 264                   TEMP = Y( K )
 265                END IF
 266                AK = A( K )
 267                PERT = SIGN( TOL, AK )
 268    70          CONTINUE
 269                ABSAK = ABS( AK )
 270                IF( ABSAK.LT.ONE ) THEN
 271                   IF( ABSAK.LT.SFMIN ) THEN
 272                      IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
 273      $                    THEN
 274                         AK = AK + PERT
 275                         PERT = 2*PERT
 276                         GO TO 70
 277                      ELSE
 278                         TEMP = TEMP*BIGNUM
 279                         AK = AK*BIGNUM
 280                      END IF
 281                   ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
 282                      AK = AK + PERT
 283                      PERT = 2*PERT
 284                      GO TO 70
 285                   END IF
 286                END IF
 287                Y( K ) = TEMP / AK
 288    80       CONTINUE
 289          END IF
 290 *
 291          DO 90 K = N, 2, -1
 292             IF( IN( K-1 ).EQ.0 ) THEN
 293                Y( K-1 ) = Y( K-1 ) - C( K-1 )*Y( K )
 294             ELSE
 295                TEMP = Y( K-1 )
 296                Y( K-1 ) = Y( K )
 297                Y( K ) = TEMP - C( K-1 )*Y( K )
 298             END IF
 299    90    CONTINUE
 300       END IF
 301 *
 302 *     End of DLAGTS
 303 *
 304       END