Return to dppcon.f CVS log

Up to [local] / rpl / lapack / lapack

Commit initial.

    1:       SUBROUTINE DPPCON( UPLO, N, AP, ANORM, RCOND, WORK, IWORK, INFO )
    2: *
    3: *  -- LAPACK routine (version 3.2) --
    4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    6: *     November 2006
    7: *
    8: *     Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH.
    9: *
   10: *     .. Scalar Arguments ..
   11:       CHARACTER          UPLO
   12:       INTEGER            INFO, N
   13:       DOUBLE PRECISION   ANORM, RCOND
   14: *     ..
   15: *     .. Array Arguments ..
   16:       INTEGER            IWORK( * )
   17:       DOUBLE PRECISION   AP( * ), WORK( * )
   18: *     ..
   19: *
   20: *  Purpose
   21: *  =======
   22: *
   23: *  DPPCON estimates the reciprocal of the condition number (in the
   24: *  1-norm) of a real symmetric positive definite packed matrix using
   25: *  the Cholesky factorization A = U**T*U or A = L*L**T computed by
   26: *  DPPTRF.
   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: *  UPLO    (input) CHARACTER*1
   35: *          = 'U':  Upper triangle of A is stored;
   36: *          = 'L':  Lower triangle of A is stored.
   37: *
   38: *  N       (input) INTEGER
   39: *          The order of the matrix A.  N >= 0.
   40: *
   41: *  AP      (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
   42: *          The triangular factor U or L from the Cholesky factorization
   43: *          A = U**T*U or A = L*L**T, packed columnwise in a linear
   44: *          array.  The j-th column of U or L is stored in the array AP
   45: *          as follows:
   46: *          if UPLO = 'U', AP(i + (j-1)*j/2) = U(i,j) for 1<=i<=j;
   47: *          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = L(i,j) for j<=i<=n.
   48: *
   49: *  ANORM   (input) DOUBLE PRECISION
   50: *          The 1-norm (or infinity-norm) of the symmetric matrix A.
   51: *
   52: *  RCOND   (output) DOUBLE PRECISION
   53: *          The reciprocal of the condition number of the matrix A,
   54: *          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
   55: *          estimate of the 1-norm of inv(A) computed in this routine.
   56: *
   57: *  WORK    (workspace) DOUBLE PRECISION array, dimension (3*N)
   58: *
   59: *  IWORK   (workspace) INTEGER array, dimension (N)
   60: *
   61: *  INFO    (output) INTEGER
   62: *          = 0:  successful exit
   63: *          < 0:  if INFO = -i, the i-th argument had an illegal value
   64: *
   65: *  =====================================================================
   66: *
   67: *     .. Parameters ..
   68:       DOUBLE PRECISION   ONE, ZERO
   69:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
   70: *     ..
   71: *     .. Local Scalars ..
   72:       LOGICAL            UPPER
   73:       CHARACTER          NORMIN
   74:       INTEGER            IX, KASE
   75:       DOUBLE PRECISION   AINVNM, SCALE, SCALEL, SCALEU, SMLNUM
   76: *     ..
   77: *     .. Local Arrays ..
   78:       INTEGER            ISAVE( 3 )
   79: *     ..
   80: *     .. External Functions ..
   81:       LOGICAL            LSAME
   82:       INTEGER            IDAMAX
   83:       DOUBLE PRECISION   DLAMCH
   84:       EXTERNAL           LSAME, IDAMAX, DLAMCH
   85: *     ..
   86: *     .. External Subroutines ..
   87:       EXTERNAL           DLACN2, DLATPS, DRSCL, XERBLA
   88: *     ..
   89: *     .. Intrinsic Functions ..
   90:       INTRINSIC          ABS
   91: *     ..
   92: *     .. Executable Statements ..
   93: *
   94: *     Test the input parameters.
   95: *
   96:       INFO = 0
   97:       UPPER = LSAME( UPLO, 'U' )
   98:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
   99:          INFO = -1
  100:       ELSE IF( N.LT.0 ) THEN
  101:          INFO = -2
  102:       ELSE IF( ANORM.LT.ZERO ) THEN
  103:          INFO = -4
  104:       END IF
  105:       IF( INFO.NE.0 ) THEN
  106:          CALL XERBLA( 'DPPCON', -INFO )
  107:          RETURN
  108:       END IF
  109: *
  110: *     Quick return if possible
  111: *
  112:       RCOND = ZERO
  113:       IF( N.EQ.0 ) THEN
  114:          RCOND = ONE
  115:          RETURN
  116:       ELSE IF( ANORM.EQ.ZERO ) THEN
  117:          RETURN
  118:       END IF
  119: *
  120:       SMLNUM = DLAMCH( 'Safe minimum' )
  121: *
  122: *     Estimate the 1-norm of the inverse.
  123: *
  124:       KASE = 0
  125:       NORMIN = 'N'
  126:    10 CONTINUE
  127:       CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
  128:       IF( KASE.NE.0 ) THEN
  129:          IF( UPPER ) THEN
  130: *
  131: *           Multiply by inv(U').
  132: *
  133:             CALL DLATPS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N,
  134:      $                   AP, WORK, SCALEL, WORK( 2*N+1 ), INFO )
  135:             NORMIN = 'Y'
  136: *
  137: *           Multiply by inv(U).
  138: *
  139:             CALL DLATPS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
  140:      $                   AP, WORK, SCALEU, WORK( 2*N+1 ), INFO )
  141:          ELSE
  142: *
  143: *           Multiply by inv(L).
  144: *
  145:             CALL DLATPS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N,
  146:      $                   AP, WORK, SCALEL, WORK( 2*N+1 ), INFO )
  147:             NORMIN = 'Y'
  148: *
  149: *           Multiply by inv(L').
  150: *
  151:             CALL DLATPS( 'Lower', 'Transpose', 'Non-unit', NORMIN, N,
  152:      $                   AP, WORK, SCALEU, WORK( 2*N+1 ), INFO )
  153:          END IF
  154: *
  155: *        Multiply by 1/SCALE if doing so will not cause overflow.
  156: *
  157:          SCALE = SCALEL*SCALEU
  158:          IF( SCALE.NE.ONE ) THEN
  159:             IX = IDAMAX( N, WORK, 1 )
  160:             IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
  161:      $         GO TO 20
  162:             CALL DRSCL( N, SCALE, WORK, 1 )
  163:          END IF
  164:          GO TO 10
  165:       END IF
  166: *
  167: *     Compute the estimate of the reciprocal condition number.
  168: *
  169:       IF( AINVNM.NE.ZERO )
  170:      $   RCOND = ( ONE / AINVNM ) / ANORM
  171: *
  172:    20 CONTINUE
  173:       RETURN
  174: *
  175: *     End of DPPCON
  176: *
  177:       END