Return to zpocon.f CVS log

Up to [local] / rpl / lapack / lapack

Commit initial.

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