Return to dgtcon.f CVS log

Up to [local] / rpl / lapack / lapack

Patches pour OS/2

    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