SRC/dlagts.f

   1 *> \brief \b 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.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DLAGTS + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlagts.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlagts.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlagts.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLAGTS( JOB, N, A, B, C, D, IN, Y, TOL, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       INTEGER            INFO, JOB, N
  25 *       DOUBLE PRECISION   TOL
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       INTEGER            IN( * )
  29 *       DOUBLE PRECISION   A( * ), B( * ), C( * ), D( * ), Y( * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> DLAGTS may be used to solve one of the systems of equations
  39 *>
  40 *>    (T - lambda*I)*x = y   or   (T - lambda*I)**T*x = y,
  41 *>
  42 *> where T is an n by n tridiagonal matrix, for x, following the
  43 *> factorization of (T - lambda*I) as
  44 *>
  45 *>    (T - lambda*I) = P*L*U ,
  46 *>
  47 *> by routine DLAGTF. The choice of equation to be solved is
  48 *> controlled by the argument JOB, and in each case there is an option
  49 *> to perturb zero or very small diagonal elements of U, this option
  50 *> being intended for use in applications such as inverse iteration.
  51 *> \endverbatim
  52 *
  53 *  Arguments:
  54 *  ==========
  55 *
  56 *> \param[in] JOB
  57 *> \verbatim
  58 *>          JOB is INTEGER
  59 *>          Specifies the job to be performed by DLAGTS as follows:
  60 *>          =  1: The equations  (T - lambda*I)x = y  are to be solved,
  61 *>                but diagonal elements of U are not to be perturbed.
  62 *>          = -1: The equations  (T - lambda*I)x = y  are to be solved
  63 *>                and, if overflow would otherwise occur, the diagonal
  64 *>                elements of U are to be perturbed. See argument TOL
  65 *>                below.
  66 *>          =  2: The equations  (T - lambda*I)**Tx = y  are to be solved,
  67 *>                but diagonal elements of U are not to be perturbed.
  68 *>          = -2: The equations  (T - lambda*I)**Tx = y  are to be solved
  69 *>                and, if overflow would otherwise occur, the diagonal
  70 *>                elements of U are to be perturbed. See argument TOL
  71 *>                below.
  72 *> \endverbatim
  73 *>
  74 *> \param[in] N
  75 *> \verbatim
  76 *>          N is INTEGER
  77 *>          The order of the matrix T.
  78 *> \endverbatim
  79 *>
  80 *> \param[in] A
  81 *> \verbatim
  82 *>          A is DOUBLE PRECISION array, dimension (N)
  83 *>          On entry, A must contain the diagonal elements of U as
  84 *>          returned from DLAGTF.
  85 *> \endverbatim
  86 *>
  87 *> \param[in] B
  88 *> \verbatim
  89 *>          B is DOUBLE PRECISION array, dimension (N-1)
  90 *>          On entry, B must contain the first super-diagonal elements of
  91 *>          U as returned from DLAGTF.
  92 *> \endverbatim
  93 *>
  94 *> \param[in] C
  95 *> \verbatim
  96 *>          C is DOUBLE PRECISION array, dimension (N-1)
  97 *>          On entry, C must contain the sub-diagonal elements of L as
  98 *>          returned from DLAGTF.
  99 *> \endverbatim
 100 *>
 101 *> \param[in] D
 102 *> \verbatim
 103 *>          D is DOUBLE PRECISION array, dimension (N-2)
 104 *>          On entry, D must contain the second super-diagonal elements
 105 *>          of U as returned from DLAGTF.
 106 *> \endverbatim
 107 *>
 108 *> \param[in] IN
 109 *> \verbatim
 110 *>          IN is INTEGER array, dimension (N)
 111 *>          On entry, IN must contain details of the matrix P as returned
 112 *>          from DLAGTF.
 113 *> \endverbatim
 114 *>
 115 *> \param[in,out] Y
 116 *> \verbatim
 117 *>          Y is DOUBLE PRECISION array, dimension (N)
 118 *>          On entry, the right hand side vector y.
 119 *>          On exit, Y is overwritten by the solution vector x.
 120 *> \endverbatim
 121 *>
 122 *> \param[in,out] TOL
 123 *> \verbatim
 124 *>          TOL is DOUBLE PRECISION
 125 *>          On entry, with  JOB .lt. 0, TOL should be the minimum
 126 *>          perturbation to be made to very small diagonal elements of U.
 127 *>          TOL should normally be chosen as about eps*norm(U), where eps
 128 *>          is the relative machine precision, but if TOL is supplied as
 129 *>          non-positive, then it is reset to eps*max( abs( u(i,j) ) ).
 130 *>          If  JOB .gt. 0  then TOL is not referenced.
 131 *>
 132 *>          On exit, TOL is changed as described above, only if TOL is
 133 *>          non-positive on entry. Otherwise TOL is unchanged.
 134 *> \endverbatim
 135 *>
 136 *> \param[out] INFO
 137 *> \verbatim
 138 *>          INFO is INTEGER
 139 *>          = 0   : successful exit
 140 *>          .lt. 0: if INFO = -i, the i-th argument had an illegal value
 141 *>          .gt. 0: overflow would occur when computing the INFO(th)
 142 *>                  element of the solution vector x. This can only occur
 143 *>                  when JOB is supplied as positive and either means
 144 *>                  that a diagonal element of U is very small, or that
 145 *>                  the elements of the right-hand side vector y are very
 146 *>                  large.
 147 *> \endverbatim
 148 *
 149 *  Authors:
 150 *  ========
 151 *
 152 *> \author Univ. of Tennessee
 153 *> \author Univ. of California Berkeley
 154 *> \author Univ. of Colorado Denver
 155 *> \author NAG Ltd.
 156 *
 157 *> \date September 2012
 158 *
 159 *> \ingroup auxOTHERauxiliary
 160 *
 161 *  =====================================================================
 162       SUBROUTINE DLAGTS( JOB, N, A, B, C, D, IN, Y, TOL, INFO )
 163 *
 164 *  -- LAPACK auxiliary routine (version 3.4.2) --
 165 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 166 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 167 *     September 2012
 168 *
 169 *     .. Scalar Arguments ..
 170       INTEGER            INFO, JOB, N
 171       DOUBLE PRECISION   TOL
 172 *     ..
 173 *     .. Array Arguments ..
 174       INTEGER            IN( * )
 175       DOUBLE PRECISION   A( * ), B( * ), C( * ), D( * ), Y( * )
 176 *     ..
 177 *
 178 *  =====================================================================
 179 *
 180 *     .. Parameters ..
 181       DOUBLE PRECISION   ONE, ZERO
 182       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 183 *     ..
 184 *     .. Local Scalars ..
 185       INTEGER            K
 186       DOUBLE PRECISION   ABSAK, AK, BIGNUM, EPS, PERT, SFMIN, TEMP
 187 *     ..
 188 *     .. Intrinsic Functions ..
 189       INTRINSIC          ABS, MAX, SIGN
 190 *     ..
 191 *     .. External Functions ..
 192       DOUBLE PRECISION   DLAMCH
 193       EXTERNAL           DLAMCH
 194 *     ..
 195 *     .. External Subroutines ..
 196       EXTERNAL           XERBLA
 197 *     ..
 198 *     .. Executable Statements ..
 199 *
 200       INFO = 0
 201       IF( ( ABS( JOB ).GT.2 ) .OR. ( JOB.EQ.0 ) ) THEN
 202          INFO = -1
 203       ELSE IF( N.LT.0 ) THEN
 204          INFO = -2
 205       END IF
 206       IF( INFO.NE.0 ) THEN
 207          CALL XERBLA( 'DLAGTS', -INFO )
 208          RETURN
 209       END IF
 210 *
 211       IF( N.EQ.0 )
 212      $   RETURN
 213 *
 214       EPS = DLAMCH( 'Epsilon' )
 215       SFMIN = DLAMCH( 'Safe minimum' )
 216       BIGNUM = ONE / SFMIN
 217 *
 218       IF( JOB.LT.0 ) THEN
 219          IF( TOL.LE.ZERO ) THEN
 220             TOL = ABS( A( 1 ) )
 221             IF( N.GT.1 )
 222      $         TOL = MAX( TOL, ABS( A( 2 ) ), ABS( B( 1 ) ) )
 223             DO 10 K = 3, N
 224                TOL = MAX( TOL, ABS( A( K ) ), ABS( B( K-1 ) ),
 225      $               ABS( D( K-2 ) ) )
 226    10       CONTINUE
 227             TOL = TOL*EPS
 228             IF( TOL.EQ.ZERO )
 229      $         TOL = EPS
 230          END IF
 231       END IF
 232 *
 233       IF( ABS( JOB ).EQ.1 ) THEN
 234          DO 20 K = 2, N
 235             IF( IN( K-1 ).EQ.0 ) THEN
 236                Y( K ) = Y( K ) - C( K-1 )*Y( K-1 )
 237             ELSE
 238                TEMP = Y( K-1 )
 239                Y( K-1 ) = Y( K )
 240                Y( K ) = TEMP - C( K-1 )*Y( K )
 241             END IF
 242    20    CONTINUE
 243          IF( JOB.EQ.1 ) THEN
 244             DO 30 K = N, 1, -1
 245                IF( K.LE.N-2 ) THEN
 246                   TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 )
 247                ELSE IF( K.EQ.N-1 ) THEN
 248                   TEMP = Y( K ) - B( K )*Y( K+1 )
 249                ELSE
 250                   TEMP = Y( K )
 251                END IF
 252                AK = A( K )
 253                ABSAK = ABS( AK )
 254                IF( ABSAK.LT.ONE ) THEN
 255                   IF( ABSAK.LT.SFMIN ) THEN
 256                      IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
 257      $                    THEN
 258                         INFO = K
 259                         RETURN
 260                      ELSE
 261                         TEMP = TEMP*BIGNUM
 262                         AK = AK*BIGNUM
 263                      END IF
 264                   ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
 265                      INFO = K
 266                      RETURN
 267                   END IF
 268                END IF
 269                Y( K ) = TEMP / AK
 270    30       CONTINUE
 271          ELSE
 272             DO 50 K = N, 1, -1
 273                IF( K.LE.N-2 ) THEN
 274                   TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 )
 275                ELSE IF( K.EQ.N-1 ) THEN
 276                   TEMP = Y( K ) - B( K )*Y( K+1 )
 277                ELSE
 278                   TEMP = Y( K )
 279                END IF
 280                AK = A( K )
 281                PERT = SIGN( TOL, AK )
 282    40          CONTINUE
 283                ABSAK = ABS( AK )
 284                IF( ABSAK.LT.ONE ) THEN
 285                   IF( ABSAK.LT.SFMIN ) THEN
 286                      IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
 287      $                    THEN
 288                         AK = AK + PERT
 289                         PERT = 2*PERT
 290                         GO TO 40
 291                      ELSE
 292                         TEMP = TEMP*BIGNUM
 293                         AK = AK*BIGNUM
 294                      END IF
 295                   ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
 296                      AK = AK + PERT
 297                      PERT = 2*PERT
 298                      GO TO 40
 299                   END IF
 300                END IF
 301                Y( K ) = TEMP / AK
 302    50       CONTINUE
 303          END IF
 304       ELSE
 305 *
 306 *        Come to here if  JOB = 2 or -2
 307 *
 308          IF( JOB.EQ.2 ) THEN
 309             DO 60 K = 1, N
 310                IF( K.GE.3 ) THEN
 311                   TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 )
 312                ELSE IF( K.EQ.2 ) THEN
 313                   TEMP = Y( K ) - B( K-1 )*Y( K-1 )
 314                ELSE
 315                   TEMP = Y( K )
 316                END IF
 317                AK = A( K )
 318                ABSAK = ABS( AK )
 319                IF( ABSAK.LT.ONE ) THEN
 320                   IF( ABSAK.LT.SFMIN ) THEN
 321                      IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
 322      $                    THEN
 323                         INFO = K
 324                         RETURN
 325                      ELSE
 326                         TEMP = TEMP*BIGNUM
 327                         AK = AK*BIGNUM
 328                      END IF
 329                   ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
 330                      INFO = K
 331                      RETURN
 332                   END IF
 333                END IF
 334                Y( K ) = TEMP / AK
 335    60       CONTINUE
 336          ELSE
 337             DO 80 K = 1, N
 338                IF( K.GE.3 ) THEN
 339                   TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 )
 340                ELSE IF( K.EQ.2 ) THEN
 341                   TEMP = Y( K ) - B( K-1 )*Y( K-1 )
 342                ELSE
 343                   TEMP = Y( K )
 344                END IF
 345                AK = A( K )
 346                PERT = SIGN( TOL, AK )
 347    70          CONTINUE
 348                ABSAK = ABS( AK )
 349                IF( ABSAK.LT.ONE ) THEN
 350                   IF( ABSAK.LT.SFMIN ) THEN
 351                      IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
 352      $                    THEN
 353                         AK = AK + PERT
 354                         PERT = 2*PERT
 355                         GO TO 70
 356                      ELSE
 357                         TEMP = TEMP*BIGNUM
 358                         AK = AK*BIGNUM
 359                      END IF
 360                   ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
 361                      AK = AK + PERT
 362                      PERT = 2*PERT
 363                      GO TO 70
 364                   END IF
 365                END IF
 366                Y( K ) = TEMP / AK
 367    80       CONTINUE
 368          END IF
 369 *
 370          DO 90 K = N, 2, -1
 371             IF( IN( K-1 ).EQ.0 ) THEN
 372                Y( K-1 ) = Y( K-1 ) - C( K-1 )*Y( K )
 373             ELSE
 374                TEMP = Y( K-1 )
 375                Y( K-1 ) = Y( K )
 376                Y( K ) = TEMP - C( K-1 )*Y( K )
 377             END IF
 378    90    CONTINUE
 379       END IF
 380 *
 381 *     End of DLAGTS
 382 *
 383       END