Annotation of rpl/lapack/lapack/dlagtf.f, revision 1.19

1.11      bertrand    1: *> \brief \b DLAGTF 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.
1.8       bertrand    2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
1.15      bertrand    5: * Online html documentation available at
                      6: *            http://www.netlib.org/lapack/explore-html/
1.8       bertrand    7: *
                      8: *> \htmlonly
1.15      bertrand    9: *> Download DLAGTF + dependencies
                     10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlagtf.f">
                     11: *> [TGZ]</a>
                     12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlagtf.f">
                     13: *> [ZIP]</a>
                     14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlagtf.f">
1.8       bertrand   15: *> [TXT]</a>
1.15      bertrand   16: *> \endhtmlonly
1.8       bertrand   17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE DLAGTF( N, A, LAMBDA, B, C, TOL, D, IN, INFO )
1.15      bertrand   22: *
1.8       bertrand   23: *       .. Scalar Arguments ..
                     24: *       INTEGER            INFO, N
                     25: *       DOUBLE PRECISION   LAMBDA, TOL
                     26: *       ..
                     27: *       .. Array Arguments ..
                     28: *       INTEGER            IN( * )
                     29: *       DOUBLE PRECISION   A( * ), B( * ), C( * ), D( * )
                     30: *       ..
1.15      bertrand   31: *
1.8       bertrand   32: *
                     33: *> \par Purpose:
                     34: *  =============
                     35: *>
                     36: *> \verbatim
                     37: *>
                     38: *> DLAGTF 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 DLAGTF may
                     52: *> be used, in conjunction with DLAGTS, 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 DOUBLE PRECISION 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 DOUBLE PRECISION
                     77: *>          On entry, the scalar lambda.
                     78: *> \endverbatim
                     79: *>
                     80: *> \param[in,out] B
                     81: *> \verbatim
                     82: *>          B is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
                    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 DOUBLE PRECISION 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: *>
1.18      bertrand  128: *>             abs( u(j,j) ) <= norm( (T - lambda*I)(j) )*TOL,
1.8       bertrand  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
1.18      bertrand  140: *>          = 0:  successful exit
                    141: *>          < 0:  if INFO = -k, the kth argument had an illegal value
1.8       bertrand  142: *> \endverbatim
                    143: *
                    144: *  Authors:
                    145: *  ========
                    146: *
1.15      bertrand  147: *> \author Univ. of Tennessee
                    148: *> \author Univ. of California Berkeley
                    149: *> \author Univ. of Colorado Denver
                    150: *> \author NAG Ltd.
1.8       bertrand  151: *
                    152: *> \ingroup auxOTHERcomputational
                    153: *
                    154: *  =====================================================================
1.1       bertrand  155:       SUBROUTINE DLAGTF( N, A, LAMBDA, B, C, TOL, D, IN, INFO )
                    156: *
1.19    ! bertrand  157: *  -- LAPACK computational routine --
1.1       bertrand  158: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    159: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                    160: *
                    161: *     .. Scalar Arguments ..
                    162:       INTEGER            INFO, N
                    163:       DOUBLE PRECISION   LAMBDA, TOL
                    164: *     ..
                    165: *     .. Array Arguments ..
                    166:       INTEGER            IN( * )
                    167:       DOUBLE PRECISION   A( * ), B( * ), C( * ), D( * )
                    168: *     ..
                    169: *
                    170: * =====================================================================
                    171: *
                    172: *     .. Parameters ..
                    173:       DOUBLE PRECISION   ZERO
                    174:       PARAMETER          ( ZERO = 0.0D+0 )
                    175: *     ..
                    176: *     .. Local Scalars ..
                    177:       INTEGER            K
                    178:       DOUBLE PRECISION   EPS, MULT, PIV1, PIV2, SCALE1, SCALE2, TEMP, TL
                    179: *     ..
                    180: *     .. Intrinsic Functions ..
                    181:       INTRINSIC          ABS, MAX
                    182: *     ..
                    183: *     .. External Functions ..
                    184:       DOUBLE PRECISION   DLAMCH
                    185:       EXTERNAL           DLAMCH
                    186: *     ..
                    187: *     .. External Subroutines ..
                    188:       EXTERNAL           XERBLA
                    189: *     ..
                    190: *     .. Executable Statements ..
                    191: *
                    192:       INFO = 0
                    193:       IF( N.LT.0 ) THEN
                    194:          INFO = -1
                    195:          CALL XERBLA( 'DLAGTF', -INFO )
                    196:          RETURN
                    197:       END IF
                    198: *
                    199:       IF( N.EQ.0 )
                    200:      $   RETURN
                    201: *
                    202:       A( 1 ) = A( 1 ) - LAMBDA
                    203:       IN( N ) = 0
                    204:       IF( N.EQ.1 ) THEN
                    205:          IF( A( 1 ).EQ.ZERO )
                    206:      $      IN( 1 ) = 1
                    207:          RETURN
                    208:       END IF
                    209: *
                    210:       EPS = DLAMCH( 'Epsilon' )
                    211: *
                    212:       TL = MAX( TOL, EPS )
                    213:       SCALE1 = ABS( A( 1 ) ) + ABS( B( 1 ) )
                    214:       DO 10 K = 1, N - 1
                    215:          A( K+1 ) = A( K+1 ) - LAMBDA
                    216:          SCALE2 = ABS( C( K ) ) + ABS( A( K+1 ) )
                    217:          IF( K.LT.( N-1 ) )
                    218:      $      SCALE2 = SCALE2 + ABS( B( K+1 ) )
                    219:          IF( A( K ).EQ.ZERO ) THEN
                    220:             PIV1 = ZERO
                    221:          ELSE
                    222:             PIV1 = ABS( A( K ) ) / SCALE1
                    223:          END IF
                    224:          IF( C( K ).EQ.ZERO ) THEN
                    225:             IN( K ) = 0
                    226:             PIV2 = ZERO
                    227:             SCALE1 = SCALE2
                    228:             IF( K.LT.( N-1 ) )
                    229:      $         D( K ) = ZERO
                    230:          ELSE
                    231:             PIV2 = ABS( C( K ) ) / SCALE2
                    232:             IF( PIV2.LE.PIV1 ) THEN
                    233:                IN( K ) = 0
                    234:                SCALE1 = SCALE2
                    235:                C( K ) = C( K ) / A( K )
                    236:                A( K+1 ) = A( K+1 ) - C( K )*B( K )
                    237:                IF( K.LT.( N-1 ) )
                    238:      $            D( K ) = ZERO
                    239:             ELSE
                    240:                IN( K ) = 1
                    241:                MULT = A( K ) / C( K )
                    242:                A( K ) = C( K )
                    243:                TEMP = A( K+1 )
                    244:                A( K+1 ) = B( K ) - MULT*TEMP
                    245:                IF( K.LT.( N-1 ) ) THEN
                    246:                   D( K ) = B( K+1 )
                    247:                   B( K+1 ) = -MULT*D( K )
                    248:                END IF
                    249:                B( K ) = TEMP
                    250:                C( K ) = MULT
                    251:             END IF
                    252:          END IF
                    253:          IF( ( MAX( PIV1, PIV2 ).LE.TL ) .AND. ( IN( N ).EQ.0 ) )
                    254:      $      IN( N ) = K
                    255:    10 CONTINUE
                    256:       IF( ( ABS( A( N ) ).LE.SCALE1*TL ) .AND. ( IN( N ).EQ.0 ) )
                    257:      $   IN( N ) = N
                    258: *
                    259:       RETURN
                    260: *
                    261: *     End of DLAGTF
                    262: *
                    263:       END

CVSweb interface <joel.bertrand@systella.fr>