Annotation of rpl/lapack/lapack/dgtcon.f, revision 1.6

1.1       bertrand    1:       SUBROUTINE DGTCON( NORM, N, DL, D, DU, DU2, IPIV, ANORM, RCOND,
                      2:      $                   WORK, IWORK, INFO )
                      3: *
                      4: *  -- LAPACK routine (version 3.2) --
                      5: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                      6: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                      7: *     November 2006
                      8: *
                      9: *     Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH.
                     10: *
                     11: *     .. Scalar Arguments ..
                     12:       CHARACTER          NORM
                     13:       INTEGER            INFO, N
                     14:       DOUBLE PRECISION   ANORM, RCOND
                     15: *     ..
                     16: *     .. Array Arguments ..
                     17:       INTEGER            IPIV( * ), IWORK( * )
                     18:       DOUBLE PRECISION   D( * ), DL( * ), DU( * ), DU2( * ), WORK( * )
                     19: *     ..
                     20: *
                     21: *  Purpose
                     22: *  =======
                     23: *
                     24: *  DGTCON estimates the reciprocal of the condition number of a real
                     25: *  tridiagonal matrix A using the LU factorization as computed by
                     26: *  DGTTRF.
                     27: *
                     28: *  An estimate is obtained for norm(inv(A)), and the reciprocal of the
                     29: *  condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
                     30: *
                     31: *  Arguments
                     32: *  =========
                     33: *
                     34: *  NORM    (input) CHARACTER*1
                     35: *          Specifies whether the 1-norm condition number or the
                     36: *          infinity-norm condition number is required:
                     37: *          = '1' or 'O':  1-norm;
                     38: *          = 'I':         Infinity-norm.
                     39: *
                     40: *  N       (input) INTEGER
                     41: *          The order of the matrix A.  N >= 0.
                     42: *
                     43: *  DL      (input) DOUBLE PRECISION array, dimension (N-1)
                     44: *          The (n-1) multipliers that define the matrix L from the
                     45: *          LU factorization of A as computed by DGTTRF.
                     46: *
                     47: *  D       (input) DOUBLE PRECISION array, dimension (N)
                     48: *          The n diagonal elements of the upper triangular matrix U from
                     49: *          the LU factorization of A.
                     50: *
                     51: *  DU      (input) DOUBLE PRECISION array, dimension (N-1)
                     52: *          The (n-1) elements of the first superdiagonal of U.
                     53: *
                     54: *  DU2     (input) DOUBLE PRECISION array, dimension (N-2)
                     55: *          The (n-2) elements of the second superdiagonal of U.
                     56: *
                     57: *  IPIV    (input) INTEGER array, dimension (N)
                     58: *          The pivot indices; for 1 <= i <= n, row i of the matrix was
                     59: *          interchanged with row IPIV(i).  IPIV(i) will always be either
                     60: *          i or i+1; IPIV(i) = i indicates a row interchange was not
                     61: *          required.
                     62: *
                     63: *  ANORM   (input) DOUBLE PRECISION
                     64: *          If NORM = '1' or 'O', the 1-norm of the original matrix A.
                     65: *          If NORM = 'I', the infinity-norm of the original matrix A.
                     66: *
                     67: *  RCOND   (output) DOUBLE PRECISION
                     68: *          The reciprocal of the condition number of the matrix A,
                     69: *          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
                     70: *          estimate of the 1-norm of inv(A) computed in this routine.
                     71: *
                     72: *  WORK    (workspace) DOUBLE PRECISION array, dimension (2*N)
                     73: *
                     74: *  IWORK   (workspace) INTEGER array, dimension (N)
                     75: *
                     76: *  INFO    (output) INTEGER
                     77: *          = 0:  successful exit
                     78: *          < 0:  if INFO = -i, the i-th argument had an illegal value
                     79: *
                     80: *  =====================================================================
                     81: *
                     82: *     .. Parameters ..
                     83:       DOUBLE PRECISION   ONE, ZERO
                     84:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
                     85: *     ..
                     86: *     .. Local Scalars ..
                     87:       LOGICAL            ONENRM
                     88:       INTEGER            I, KASE, KASE1
                     89:       DOUBLE PRECISION   AINVNM
                     90: *     ..
                     91: *     .. Local Arrays ..
                     92:       INTEGER            ISAVE( 3 )
                     93: *     ..
                     94: *     .. External Functions ..
                     95:       LOGICAL            LSAME
                     96:       EXTERNAL           LSAME
                     97: *     ..
                     98: *     .. External Subroutines ..
                     99:       EXTERNAL           DGTTRS, DLACN2, XERBLA
                    100: *     ..
                    101: *     .. Executable Statements ..
                    102: *
                    103: *     Test the input arguments.
                    104: *
                    105:       INFO = 0
                    106:       ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
                    107:       IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
                    108:          INFO = -1
                    109:       ELSE IF( N.LT.0 ) THEN
                    110:          INFO = -2
                    111:       ELSE IF( ANORM.LT.ZERO ) THEN
                    112:          INFO = -8
                    113:       END IF
                    114:       IF( INFO.NE.0 ) THEN
                    115:          CALL XERBLA( 'DGTCON', -INFO )
                    116:          RETURN
                    117:       END IF
                    118: *
                    119: *     Quick return if possible
                    120: *
                    121:       RCOND = ZERO
                    122:       IF( N.EQ.0 ) THEN
                    123:          RCOND = ONE
                    124:          RETURN
                    125:       ELSE IF( ANORM.EQ.ZERO ) THEN
                    126:          RETURN
                    127:       END IF
                    128: *
                    129: *     Check that D(1:N) is non-zero.
                    130: *
                    131:       DO 10 I = 1, N
                    132:          IF( D( I ).EQ.ZERO )
                    133:      $      RETURN
                    134:    10 CONTINUE
                    135: *
                    136:       AINVNM = ZERO
                    137:       IF( ONENRM ) THEN
                    138:          KASE1 = 1
                    139:       ELSE
                    140:          KASE1 = 2
                    141:       END IF
                    142:       KASE = 0
                    143:    20 CONTINUE
                    144:       CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
                    145:       IF( KASE.NE.0 ) THEN
                    146:          IF( KASE.EQ.KASE1 ) THEN
                    147: *
                    148: *           Multiply by inv(U)*inv(L).
                    149: *
                    150:             CALL DGTTRS( 'No transpose', N, 1, DL, D, DU, DU2, IPIV,
                    151:      $                   WORK, N, INFO )
                    152:          ELSE
                    153: *
                    154: *           Multiply by inv(L')*inv(U').
                    155: *
                    156:             CALL DGTTRS( 'Transpose', N, 1, DL, D, DU, DU2, IPIV, WORK,
                    157:      $                   N, INFO )
                    158:          END IF
                    159:          GO TO 20
                    160:       END IF
                    161: *
                    162: *     Compute the estimate of the reciprocal condition number.
                    163: *
                    164:       IF( AINVNM.NE.ZERO )
                    165:      $   RCOND = ( ONE / AINVNM ) / ANORM
                    166: *
                    167:       RETURN
                    168: *
                    169: *     End of DGTCON
                    170: *
                    171:       END

CVSweb interface <joel.bertrand@systella.fr>