Return to zppcon.f CVS log

Up to [local] / rpl / lapack / lapack

Commit initial.

    1:       SUBROUTINE ZPPCON( UPLO, N, AP, ANORM, RCOND, WORK, RWORK, 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 ZLACN2 in place of ZLACON, 10 Feb 03, SJH.
    9: *
   10: *     .. Scalar Arguments ..
   11:       CHARACTER          UPLO
   12:       INTEGER            INFO, N
   13:       DOUBLE PRECISION   ANORM, RCOND
   14: *     ..
   15: *     .. Array Arguments ..
   16:       DOUBLE PRECISION   RWORK( * )
   17:       COMPLEX*16         AP( * ), WORK( * )
   18: *     ..
   19: *
   20: *  Purpose
   21: *  =======
   22: *
   23: *  ZPPCON estimates the reciprocal of the condition number (in the
   24: *  1-norm) of a complex Hermitian positive definite packed matrix using
   25: *  the Cholesky factorization A = U**H*U or A = L*L**H computed by
   26: *  ZPPTRF.
   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) COMPLEX*16 array, dimension (N*(N+1)/2)
   42: *          The triangular factor U or L from the Cholesky factorization
   43: *          A = U**H*U or A = L*L**H, 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 Hermitian 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) COMPLEX*16 array, dimension (2*N)
   58: *
   59: *  RWORK   (workspace) DOUBLE PRECISION 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:       COMPLEX*16         ZDUM
   77: *     ..
   78: *     .. Local Arrays ..
   79:       INTEGER            ISAVE( 3 )
   80: *     ..
   81: *     .. External Functions ..
   82:       LOGICAL            LSAME
   83:       INTEGER            IZAMAX
   84:       DOUBLE PRECISION   DLAMCH
   85:       EXTERNAL           LSAME, IZAMAX, DLAMCH
   86: *     ..
   87: *     .. External Subroutines ..
   88:       EXTERNAL           XERBLA, ZDRSCL, ZLACN2, ZLATPS
   89: *     ..
   90: *     .. Intrinsic Functions ..
   91:       INTRINSIC          ABS, DBLE, DIMAG
   92: *     ..
   93: *     .. Statement Functions ..
   94:       DOUBLE PRECISION   CABS1
   95: *     ..
   96: *     .. Statement Function definitions ..
   97:       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
   98: *     ..
   99: *     .. Executable Statements ..
  100: *
  101: *     Test the input parameters.
  102: *
  103:       INFO = 0
  104:       UPPER = LSAME( UPLO, 'U' )
  105:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  106:          INFO = -1
  107:       ELSE IF( N.LT.0 ) THEN
  108:          INFO = -2
  109:       ELSE IF( ANORM.LT.ZERO ) THEN
  110:          INFO = -4
  111:       END IF
  112:       IF( INFO.NE.0 ) THEN
  113:          CALL XERBLA( 'ZPPCON', -INFO )
  114:          RETURN
  115:       END IF
  116: *
  117: *     Quick return if possible
  118: *
  119:       RCOND = ZERO
  120:       IF( N.EQ.0 ) THEN
  121:          RCOND = ONE
  122:          RETURN
  123:       ELSE IF( ANORM.EQ.ZERO ) THEN
  124:          RETURN
  125:       END IF
  126: *
  127:       SMLNUM = DLAMCH( 'Safe minimum' )
  128: *
  129: *     Estimate the 1-norm of the inverse.
  130: *
  131:       KASE = 0
  132:       NORMIN = 'N'
  133:    10 CONTINUE
  134:       CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
  135:       IF( KASE.NE.0 ) THEN
  136:          IF( UPPER ) THEN
  137: *
  138: *           Multiply by inv(U').
  139: *
  140:             CALL ZLATPS( 'Upper', 'Conjugate transpose', 'Non-unit',
  141:      $                   NORMIN, N, AP, WORK, SCALEL, RWORK, INFO )
  142:             NORMIN = 'Y'
  143: *
  144: *           Multiply by inv(U).
  145: *
  146:             CALL ZLATPS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
  147:      $                   AP, WORK, SCALEU, RWORK, INFO )
  148:          ELSE
  149: *
  150: *           Multiply by inv(L).
  151: *
  152:             CALL ZLATPS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N,
  153:      $                   AP, WORK, SCALEL, RWORK, INFO )
  154:             NORMIN = 'Y'
  155: *
  156: *           Multiply by inv(L').
  157: *
  158:             CALL ZLATPS( 'Lower', 'Conjugate transpose', 'Non-unit',
  159:      $                   NORMIN, N, AP, WORK, SCALEU, RWORK, INFO )
  160:          END IF
  161: *
  162: *        Multiply by 1/SCALE if doing so will not cause overflow.
  163: *
  164:          SCALE = SCALEL*SCALEU
  165:          IF( SCALE.NE.ONE ) THEN
  166:             IX = IZAMAX( N, WORK, 1 )
  167:             IF( SCALE.LT.CABS1( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
  168:      $         GO TO 20
  169:             CALL ZDRSCL( N, SCALE, WORK, 1 )
  170:          END IF
  171:          GO TO 10
  172:       END IF
  173: *
  174: *     Compute the estimate of the reciprocal condition number.
  175: *
  176:       IF( AINVNM.NE.ZERO )
  177:      $   RCOND = ( ONE / AINVNM ) / ANORM
  178: *
  179:    20 CONTINUE
  180:       RETURN
  181: *
  182: *     End of ZPPCON
  183: *
  184:       END