LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
slagts.f
Go to the documentation of this file.
00001 *> \brief \b SLAGTS
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SLAGTS + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slagts.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slagts.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slagts.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SLAGTS( JOB, N, A, B, C, D, IN, Y, TOL, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       INTEGER            INFO, JOB, N
00025 *       REAL               TOL
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       INTEGER            IN( * )
00029 *       REAL               A( * ), B( * ), C( * ), D( * ), Y( * )
00030 *       ..
00031 *  
00032 *
00033 *> \par Purpose:
00034 *  =============
00035 *>
00036 *> \verbatim
00037 *>
00038 *> SLAGTS may be used to solve one of the systems of equations
00039 *>
00040 *>    (T - lambda*I)*x = y   or   (T - lambda*I)**T*x = y,
00041 *>
00042 *> where T is an n by n tridiagonal matrix, for x, following the
00043 *> factorization of (T - lambda*I) as
00044 *>
00045 *>    (T - lambda*I) = P*L*U ,
00046 *>
00047 *> by routine SLAGTF. The choice of equation to be solved is
00048 *> controlled by the argument JOB, and in each case there is an option
00049 *> to perturb zero or very small diagonal elements of U, this option
00050 *> being intended for use in applications such as inverse iteration.
00051 *> \endverbatim
00052 *
00053 *  Arguments:
00054 *  ==========
00055 *
00056 *> \param[in] JOB
00057 *> \verbatim
00058 *>          JOB is INTEGER
00059 *>          Specifies the job to be performed by SLAGTS as follows:
00060 *>          =  1: The equations  (T - lambda*I)x = y  are to be solved,
00061 *>                but diagonal elements of U are not to be perturbed.
00062 *>          = -1: The equations  (T - lambda*I)x = y  are to be solved
00063 *>                and, if overflow would otherwise occur, the diagonal
00064 *>                elements of U are to be perturbed. See argument TOL
00065 *>                below.
00066 *>          =  2: The equations  (T - lambda*I)**Tx = y  are to be solved,
00067 *>                but diagonal elements of U are not to be perturbed.
00068 *>          = -2: The equations  (T - lambda*I)**Tx = y  are to be solved
00069 *>                and, if overflow would otherwise occur, the diagonal
00070 *>                elements of U are to be perturbed. See argument TOL
00071 *>                below.
00072 *> \endverbatim
00073 *>
00074 *> \param[in] N
00075 *> \verbatim
00076 *>          N is INTEGER
00077 *>          The order of the matrix T.
00078 *> \endverbatim
00079 *>
00080 *> \param[in] A
00081 *> \verbatim
00082 *>          A is REAL array, dimension (N)
00083 *>          On entry, A must contain the diagonal elements of U as
00084 *>          returned from SLAGTF.
00085 *> \endverbatim
00086 *>
00087 *> \param[in] B
00088 *> \verbatim
00089 *>          B is REAL array, dimension (N-1)
00090 *>          On entry, B must contain the first super-diagonal elements of
00091 *>          U as returned from SLAGTF.
00092 *> \endverbatim
00093 *>
00094 *> \param[in] C
00095 *> \verbatim
00096 *>          C is REAL array, dimension (N-1)
00097 *>          On entry, C must contain the sub-diagonal elements of L as
00098 *>          returned from SLAGTF.
00099 *> \endverbatim
00100 *>
00101 *> \param[in] D
00102 *> \verbatim
00103 *>          D is REAL array, dimension (N-2)
00104 *>          On entry, D must contain the second super-diagonal elements
00105 *>          of U as returned from SLAGTF.
00106 *> \endverbatim
00107 *>
00108 *> \param[in] IN
00109 *> \verbatim
00110 *>          IN is INTEGER array, dimension (N)
00111 *>          On entry, IN must contain details of the matrix P as returned
00112 *>          from SLAGTF.
00113 *> \endverbatim
00114 *>
00115 *> \param[in,out] Y
00116 *> \verbatim
00117 *>          Y is REAL array, dimension (N)
00118 *>          On entry, the right hand side vector y.
00119 *>          On exit, Y is overwritten by the solution vector x.
00120 *> \endverbatim
00121 *>
00122 *> \param[in,out] TOL
00123 *> \verbatim
00124 *>          TOL is REAL
00125 *>          On entry, with  JOB .lt. 0, TOL should be the minimum
00126 *>          perturbation to be made to very small diagonal elements of U.
00127 *>          TOL should normally be chosen as about eps*norm(U), where eps
00128 *>          is the relative machine precision, but if TOL is supplied as
00129 *>          non-positive, then it is reset to eps*max( abs( u(i,j) ) ).
00130 *>          If  JOB .gt. 0  then TOL is not referenced.
00131 *>
00132 *>          On exit, TOL is changed as described above, only if TOL is
00133 *>          non-positive on entry. Otherwise TOL is unchanged.
00134 *> \endverbatim
00135 *>
00136 *> \param[out] INFO
00137 *> \verbatim
00138 *>          INFO is INTEGER
00139 *>          = 0   : successful exit
00140 *>          .lt. 0: if INFO = -i, the i-th argument had an illegal value
00141 *>          .gt. 0: overflow would occur when computing the INFO(th)
00142 *>                  element of the solution vector x. This can only occur
00143 *>                  when JOB is supplied as positive and either means
00144 *>                  that a diagonal element of U is very small, or that
00145 *>                  the elements of the right-hand side vector y are very
00146 *>                  large.
00147 *> \endverbatim
00148 *
00149 *  Authors:
00150 *  ========
00151 *
00152 *> \author Univ. of Tennessee 
00153 *> \author Univ. of California Berkeley 
00154 *> \author Univ. of Colorado Denver 
00155 *> \author NAG Ltd. 
00156 *
00157 *> \date November 2011
00158 *
00159 *> \ingroup auxOTHERauxiliary
00160 *
00161 *  =====================================================================
00162       SUBROUTINE SLAGTS( JOB, N, A, B, C, D, IN, Y, TOL, INFO )
00163 *
00164 *  -- LAPACK auxiliary routine (version 3.4.0) --
00165 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00166 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00167 *     November 2011
00168 *
00169 *     .. Scalar Arguments ..
00170       INTEGER            INFO, JOB, N
00171       REAL               TOL
00172 *     ..
00173 *     .. Array Arguments ..
00174       INTEGER            IN( * )
00175       REAL               A( * ), B( * ), C( * ), D( * ), Y( * )
00176 *     ..
00177 *
00178 *  =====================================================================
00179 *
00180 *     .. Parameters ..
00181       REAL               ONE, ZERO
00182       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
00183 *     ..
00184 *     .. Local Scalars ..
00185       INTEGER            K
00186       REAL               ABSAK, AK, BIGNUM, EPS, PERT, SFMIN, TEMP
00187 *     ..
00188 *     .. Intrinsic Functions ..
00189       INTRINSIC          ABS, MAX, SIGN
00190 *     ..
00191 *     .. External Functions ..
00192       REAL               SLAMCH
00193       EXTERNAL           SLAMCH
00194 *     ..
00195 *     .. External Subroutines ..
00196       EXTERNAL           XERBLA
00197 *     ..
00198 *     .. Executable Statements ..
00199 *
00200       INFO = 0
00201       IF( ( ABS( JOB ).GT.2 ) .OR. ( JOB.EQ.0 ) ) THEN
00202          INFO = -1
00203       ELSE IF( N.LT.0 ) THEN
00204          INFO = -2
00205       END IF
00206       IF( INFO.NE.0 ) THEN
00207          CALL XERBLA( 'SLAGTS', -INFO )
00208          RETURN
00209       END IF
00210 *
00211       IF( N.EQ.0 )
00212      $   RETURN
00213 *
00214       EPS = SLAMCH( 'Epsilon' )
00215       SFMIN = SLAMCH( 'Safe minimum' )
00216       BIGNUM = ONE / SFMIN
00217 *
00218       IF( JOB.LT.0 ) THEN
00219          IF( TOL.LE.ZERO ) THEN
00220             TOL = ABS( A( 1 ) )
00221             IF( N.GT.1 )
00222      $         TOL = MAX( TOL, ABS( A( 2 ) ), ABS( B( 1 ) ) )
00223             DO 10 K = 3, N
00224                TOL = MAX( TOL, ABS( A( K ) ), ABS( B( K-1 ) ),
00225      $               ABS( D( K-2 ) ) )
00226    10       CONTINUE
00227             TOL = TOL*EPS
00228             IF( TOL.EQ.ZERO )
00229      $         TOL = EPS
00230          END IF
00231       END IF
00232 *
00233       IF( ABS( JOB ).EQ.1 ) THEN
00234          DO 20 K = 2, N
00235             IF( IN( K-1 ).EQ.0 ) THEN
00236                Y( K ) = Y( K ) - C( K-1 )*Y( K-1 )
00237             ELSE
00238                TEMP = Y( K-1 )
00239                Y( K-1 ) = Y( K )
00240                Y( K ) = TEMP - C( K-1 )*Y( K )
00241             END IF
00242    20    CONTINUE
00243          IF( JOB.EQ.1 ) THEN
00244             DO 30 K = N, 1, -1
00245                IF( K.LE.N-2 ) THEN
00246                   TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 )
00247                ELSE IF( K.EQ.N-1 ) THEN
00248                   TEMP = Y( K ) - B( K )*Y( K+1 )
00249                ELSE
00250                   TEMP = Y( K )
00251                END IF
00252                AK = A( K )
00253                ABSAK = ABS( AK )
00254                IF( ABSAK.LT.ONE ) THEN
00255                   IF( ABSAK.LT.SFMIN ) THEN
00256                      IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
00257      $                    THEN
00258                         INFO = K
00259                         RETURN
00260                      ELSE
00261                         TEMP = TEMP*BIGNUM
00262                         AK = AK*BIGNUM
00263                      END IF
00264                   ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
00265                      INFO = K
00266                      RETURN
00267                   END IF
00268                END IF
00269                Y( K ) = TEMP / AK
00270    30       CONTINUE
00271          ELSE
00272             DO 50 K = N, 1, -1
00273                IF( K.LE.N-2 ) THEN
00274                   TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 )
00275                ELSE IF( K.EQ.N-1 ) THEN
00276                   TEMP = Y( K ) - B( K )*Y( K+1 )
00277                ELSE
00278                   TEMP = Y( K )
00279                END IF
00280                AK = A( K )
00281                PERT = SIGN( TOL, AK )
00282    40          CONTINUE
00283                ABSAK = ABS( AK )
00284                IF( ABSAK.LT.ONE ) THEN
00285                   IF( ABSAK.LT.SFMIN ) THEN
00286                      IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
00287      $                    THEN
00288                         AK = AK + PERT
00289                         PERT = 2*PERT
00290                         GO TO 40
00291                      ELSE
00292                         TEMP = TEMP*BIGNUM
00293                         AK = AK*BIGNUM
00294                      END IF
00295                   ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
00296                      AK = AK + PERT
00297                      PERT = 2*PERT
00298                      GO TO 40
00299                   END IF
00300                END IF
00301                Y( K ) = TEMP / AK
00302    50       CONTINUE
00303          END IF
00304       ELSE
00305 *
00306 *        Come to here if  JOB = 2 or -2
00307 *
00308          IF( JOB.EQ.2 ) THEN
00309             DO 60 K = 1, N
00310                IF( K.GE.3 ) THEN
00311                   TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 )
00312                ELSE IF( K.EQ.2 ) THEN
00313                   TEMP = Y( K ) - B( K-1 )*Y( K-1 )
00314                ELSE
00315                   TEMP = Y( K )
00316                END IF
00317                AK = A( K )
00318                ABSAK = ABS( AK )
00319                IF( ABSAK.LT.ONE ) THEN
00320                   IF( ABSAK.LT.SFMIN ) THEN
00321                      IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
00322      $                    THEN
00323                         INFO = K
00324                         RETURN
00325                      ELSE
00326                         TEMP = TEMP*BIGNUM
00327                         AK = AK*BIGNUM
00328                      END IF
00329                   ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
00330                      INFO = K
00331                      RETURN
00332                   END IF
00333                END IF
00334                Y( K ) = TEMP / AK
00335    60       CONTINUE
00336          ELSE
00337             DO 80 K = 1, N
00338                IF( K.GE.3 ) THEN
00339                   TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 )
00340                ELSE IF( K.EQ.2 ) THEN
00341                   TEMP = Y( K ) - B( K-1 )*Y( K-1 )
00342                ELSE
00343                   TEMP = Y( K )
00344                END IF
00345                AK = A( K )
00346                PERT = SIGN( TOL, AK )
00347    70          CONTINUE
00348                ABSAK = ABS( AK )
00349                IF( ABSAK.LT.ONE ) THEN
00350                   IF( ABSAK.LT.SFMIN ) THEN
00351                      IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
00352      $                    THEN
00353                         AK = AK + PERT
00354                         PERT = 2*PERT
00355                         GO TO 70
00356                      ELSE
00357                         TEMP = TEMP*BIGNUM
00358                         AK = AK*BIGNUM
00359                      END IF
00360                   ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
00361                      AK = AK + PERT
00362                      PERT = 2*PERT
00363                      GO TO 70
00364                   END IF
00365                END IF
00366                Y( K ) = TEMP / AK
00367    80       CONTINUE
00368          END IF
00369 *
00370          DO 90 K = N, 2, -1
00371             IF( IN( K-1 ).EQ.0 ) THEN
00372                Y( K-1 ) = Y( K-1 ) - C( K-1 )*Y( K )
00373             ELSE
00374                TEMP = Y( K-1 )
00375                Y( K-1 ) = Y( K )
00376                Y( K ) = TEMP - C( K-1 )*Y( K )
00377             END IF
00378    90    CONTINUE
00379       END IF
00380 *
00381 *     End of SLAGTS
00382 *
00383       END
 All Files Functions