SRC/slagtf.f

   1 *> \brief \b SLAGTF computes an LU factorization of a matrix T-λI, where T is a general tridiagonal matrix, and λ a scalar, using partial pivoting with row interchanges.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SLAGTF + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slagtf.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slagtf.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slagtf.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SLAGTF( N, A, LAMBDA, B, C, TOL, D, IN, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       INTEGER            INFO, N
  25 *       REAL               LAMBDA, TOL
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       INTEGER            IN( * )
  29 *       REAL               A( * ), B( * ), C( * ), D( * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> SLAGTF factorizes the matrix (T - lambda*I), where T is an n by n
  39 *> tridiagonal matrix and lambda is a scalar, as
  40 *>
  41 *>    T - lambda*I = PLU,
  42 *>
  43 *> where P is a permutation matrix, L is a unit lower tridiagonal matrix
  44 *> with at most one non-zero sub-diagonal elements per column and U is
  45 *> an upper triangular matrix with at most two non-zero super-diagonal
  46 *> elements per column.
  47 *>
  48 *> The factorization is obtained by Gaussian elimination with partial
  49 *> pivoting and implicit row scaling.
  50 *>
  51 *> The parameter LAMBDA is included in the routine so that SLAGTF may
  52 *> be used, in conjunction with SLAGTS, to obtain eigenvectors of T by
  53 *> inverse iteration.
  54 *> \endverbatim
  55 *
  56 *  Arguments:
  57 *  ==========
  58 *
  59 *> \param[in] N
  60 *> \verbatim
  61 *>          N is INTEGER
  62 *>          The order of the matrix T.
  63 *> \endverbatim
  64 *>
  65 *> \param[in,out] A
  66 *> \verbatim
  67 *>          A is REAL array, dimension (N)
  68 *>          On entry, A must contain the diagonal elements of T.
  69 *>
  70 *>          On exit, A is overwritten by the n diagonal elements of the
  71 *>          upper triangular matrix U of the factorization of T.
  72 *> \endverbatim
  73 *>
  74 *> \param[in] LAMBDA
  75 *> \verbatim
  76 *>          LAMBDA is REAL
  77 *>          On entry, the scalar lambda.
  78 *> \endverbatim
  79 *>
  80 *> \param[in,out] B
  81 *> \verbatim
  82 *>          B is REAL array, dimension (N-1)
  83 *>          On entry, B must contain the (n-1) super-diagonal elements of
  84 *>          T.
  85 *>
  86 *>          On exit, B is overwritten by the (n-1) super-diagonal
  87 *>          elements of the matrix U of the factorization of T.
  88 *> \endverbatim
  89 *>
  90 *> \param[in,out] C
  91 *> \verbatim
  92 *>          C is REAL array, dimension (N-1)
  93 *>          On entry, C must contain the (n-1) sub-diagonal elements of
  94 *>          T.
  95 *>
  96 *>          On exit, C is overwritten by the (n-1) sub-diagonal elements
  97 *>          of the matrix L of the factorization of T.
  98 *> \endverbatim
  99 *>
 100 *> \param[in] TOL
 101 *> \verbatim
 102 *>          TOL is REAL
 103 *>          On entry, a relative tolerance used to indicate whether or
 104 *>          not the matrix (T - lambda*I) is nearly singular. TOL should
 105 *>          normally be chose as approximately the largest relative error
 106 *>          in the elements of T. For example, if the elements of T are
 107 *>          correct to about 4 significant figures, then TOL should be
 108 *>          set to about 5*10**(-4). If TOL is supplied as less than eps,
 109 *>          where eps is the relative machine precision, then the value
 110 *>          eps is used in place of TOL.
 111 *> \endverbatim
 112 *>
 113 *> \param[out] D
 114 *> \verbatim
 115 *>          D is REAL array, dimension (N-2)
 116 *>          On exit, D is overwritten by the (n-2) second super-diagonal
 117 *>          elements of the matrix U of the factorization of T.
 118 *> \endverbatim
 119 *>
 120 *> \param[out] IN
 121 *> \verbatim
 122 *>          IN is INTEGER array, dimension (N)
 123 *>          On exit, IN contains details of the permutation matrix P. If
 124 *>          an interchange occurred at the kth step of the elimination,
 125 *>          then IN(k) = 1, otherwise IN(k) = 0. The element IN(n)
 126 *>          returns the smallest positive integer j such that
 127 *>
 128 *>             abs( u(j,j) ).le. norm( (T - lambda*I)(j) )*TOL,
 129 *>
 130 *>          where norm( A(j) ) denotes the sum of the absolute values of
 131 *>          the jth row of the matrix A. If no such j exists then IN(n)
 132 *>          is returned as zero. If IN(n) is returned as positive, then a
 133 *>          diagonal element of U is small, indicating that
 134 *>          (T - lambda*I) is singular or nearly singular,
 135 *> \endverbatim
 136 *>
 137 *> \param[out] INFO
 138 *> \verbatim
 139 *>          INFO is INTEGER
 140 *>          = 0   : successful exit
 141 *>          .lt. 0: if INFO = -k, the kth argument had an illegal value
 142 *> \endverbatim
 143 *
 144 *  Authors:
 145 *  ========
 146 *
 147 *> \author Univ. of Tennessee
 148 *> \author Univ. of California Berkeley
 149 *> \author Univ. of Colorado Denver
 150 *> \author NAG Ltd.
 151 *
 152 *> \date September 2012
 153 *
 154 *> \ingroup auxOTHERcomputational
 155 *
 156 *  =====================================================================
 157       SUBROUTINE SLAGTF( N, A, LAMBDA, B, C, TOL, D, IN, INFO )
 158 *
 159 *  -- LAPACK computational routine (version 3.4.2) --
 160 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 161 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 162 *     September 2012
 163 *
 164 *     .. Scalar Arguments ..
 165       INTEGER            INFO, N
 166       REAL               LAMBDA, TOL
 167 *     ..
 168 *     .. Array Arguments ..
 169       INTEGER            IN( * )
 170       REAL               A( * ), B( * ), C( * ), D( * )
 171 *     ..
 172 *
 173 * =====================================================================
 174 *
 175 *     .. Parameters ..
 176       REAL               ZERO
 177       PARAMETER          ( ZERO = 0.0E+0 )
 178 *     ..
 179 *     .. Local Scalars ..
 180       INTEGER            K
 181       REAL               EPS, MULT, PIV1, PIV2, SCALE1, SCALE2, TEMP, TL
 182 *     ..
 183 *     .. Intrinsic Functions ..
 184       INTRINSIC          ABS, MAX
 185 *     ..
 186 *     .. External Functions ..
 187       REAL               SLAMCH
 188       EXTERNAL           SLAMCH
 189 *     ..
 190 *     .. External Subroutines ..
 191       EXTERNAL           XERBLA
 192 *     ..
 193 *     .. Executable Statements ..
 194 *
 195       INFO = 0
 196       IF( N.LT.0 ) THEN
 197          INFO = -1
 198          CALL XERBLA( 'SLAGTF', -INFO )
 199          RETURN
 200       END IF
 201 *
 202       IF( N.EQ.0 )
 203      $   RETURN
 204 *
 205       A( 1 ) = A( 1 ) - LAMBDA
 206       IN( N ) = 0
 207       IF( N.EQ.1 ) THEN
 208          IF( A( 1 ).EQ.ZERO )
 209      $      IN( 1 ) = 1
 210          RETURN
 211       END IF
 212 *
 213       EPS = SLAMCH( 'Epsilon' )
 214 *
 215       TL = MAX( TOL, EPS )
 216       SCALE1 = ABS( A( 1 ) ) + ABS( B( 1 ) )
 217       DO 10 K = 1, N - 1
 218          A( K+1 ) = A( K+1 ) - LAMBDA
 219          SCALE2 = ABS( C( K ) ) + ABS( A( K+1 ) )
 220          IF( K.LT.( N-1 ) )
 221      $      SCALE2 = SCALE2 + ABS( B( K+1 ) )
 222          IF( A( K ).EQ.ZERO ) THEN
 223             PIV1 = ZERO
 224          ELSE
 225             PIV1 = ABS( A( K ) ) / SCALE1
 226          END IF
 227          IF( C( K ).EQ.ZERO ) THEN
 228             IN( K ) = 0
 229             PIV2 = ZERO
 230             SCALE1 = SCALE2
 231             IF( K.LT.( N-1 ) )
 232      $         D( K ) = ZERO
 233          ELSE
 234             PIV2 = ABS( C( K ) ) / SCALE2
 235             IF( PIV2.LE.PIV1 ) THEN
 236                IN( K ) = 0
 237                SCALE1 = SCALE2
 238                C( K ) = C( K ) / A( K )
 239                A( K+1 ) = A( K+1 ) - C( K )*B( K )
 240                IF( K.LT.( N-1 ) )
 241      $            D( K ) = ZERO
 242             ELSE
 243                IN( K ) = 1
 244                MULT = A( K ) / C( K )
 245                A( K ) = C( K )
 246                TEMP = A( K+1 )
 247                A( K+1 ) = B( K ) - MULT*TEMP
 248                IF( K.LT.( N-1 ) ) THEN
 249                   D( K ) = B( K+1 )
 250                   B( K+1 ) = -MULT*D( K )
 251                END IF
 252                B( K ) = TEMP
 253                C( K ) = MULT
 254             END IF
 255          END IF
 256          IF( ( MAX( PIV1, PIV2 ).LE.TL ) .AND. ( IN( N ).EQ.0 ) )
 257      $      IN( N ) = K
 258    10 CONTINUE
 259       IF( ( ABS( A( N ) ).LE.SCALE1*TL ) .AND. ( IN( N ).EQ.0 ) )
 260      $   IN( N ) = N
 261 *
 262       RETURN
 263 *
 264 *     End of SLAGTF
 265 *
 266       END