Return to zspcon.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

    1:       SUBROUTINE ZSPCON( UPLO, N, AP, IPIV, ANORM, RCOND, WORK, 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:       INTEGER            IPIV( * )
   17:       COMPLEX*16         AP( * ), WORK( * )
   18: *     ..
   19: *
   20: *  Purpose
   21: *  =======
   22: *
   23: *  ZSPCON estimates the reciprocal of the condition number (in the
   24: *  1-norm) of a complex symmetric packed matrix A using the
   25: *  factorization A = U*D*U**T or A = L*D*L**T computed by ZSPTRF.
   26: *
   27: *  An estimate is obtained for norm(inv(A)), and the reciprocal of the
   28: *  condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
   29: *
   30: *  Arguments
   31: *  =========
   32: *
   33: *  UPLO    (input) CHARACTER*1
   34: *          Specifies whether the details of the factorization are stored
   35: *          as an upper or lower triangular matrix.
   36: *          = 'U':  Upper triangular, form is A = U*D*U**T;
   37: *          = 'L':  Lower triangular, form is A = L*D*L**T.
   38: *
   39: *  N       (input) INTEGER
   40: *          The order of the matrix A.  N >= 0.
   41: *
   42: *  AP      (input) COMPLEX*16 array, dimension (N*(N+1)/2)
   43: *          The block diagonal matrix D and the multipliers used to
   44: *          obtain the factor U or L as computed by ZSPTRF, stored as a
   45: *          packed triangular matrix.
   46: *
   47: *  IPIV    (input) INTEGER array, dimension (N)
   48: *          Details of the interchanges and the block structure of D
   49: *          as determined by ZSPTRF.
   50: *
   51: *  ANORM   (input) DOUBLE PRECISION
   52: *          The 1-norm of the original matrix A.
   53: *
   54: *  RCOND   (output) DOUBLE PRECISION
   55: *          The reciprocal of the condition number of the matrix A,
   56: *          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
   57: *          estimate of the 1-norm of inv(A) computed in this routine.
   58: *
   59: *  WORK    (workspace) COMPLEX*16 array, dimension (2*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:       INTEGER            I, IP, KASE
   74:       DOUBLE PRECISION   AINVNM
   75: *     ..
   76: *     .. Local Arrays ..
   77:       INTEGER            ISAVE( 3 )
   78: *     ..
   79: *     .. External Functions ..
   80:       LOGICAL            LSAME
   81:       EXTERNAL           LSAME
   82: *     ..
   83: *     .. External Subroutines ..
   84:       EXTERNAL           XERBLA, ZLACN2, ZSPTRS
   85: *     ..
   86: *     .. Executable Statements ..
   87: *
   88: *     Test the input parameters.
   89: *
   90:       INFO = 0
   91:       UPPER = LSAME( UPLO, 'U' )
   92:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
   93:          INFO = -1
   94:       ELSE IF( N.LT.0 ) THEN
   95:          INFO = -2
   96:       ELSE IF( ANORM.LT.ZERO ) THEN
   97:          INFO = -5
   98:       END IF
   99:       IF( INFO.NE.0 ) THEN
  100:          CALL XERBLA( 'ZSPCON', -INFO )
  101:          RETURN
  102:       END IF
  103: *
  104: *     Quick return if possible
  105: *
  106:       RCOND = ZERO
  107:       IF( N.EQ.0 ) THEN
  108:          RCOND = ONE
  109:          RETURN
  110:       ELSE IF( ANORM.LE.ZERO ) THEN
  111:          RETURN
  112:       END IF
  113: *
  114: *     Check that the diagonal matrix D is nonsingular.
  115: *
  116:       IF( UPPER ) THEN
  117: *
  118: *        Upper triangular storage: examine D from bottom to top
  119: *
  120:          IP = N*( N+1 ) / 2
  121:          DO 10 I = N, 1, -1
  122:             IF( IPIV( I ).GT.0 .AND. AP( IP ).EQ.ZERO )
  123:      $         RETURN
  124:             IP = IP - I
  125:    10    CONTINUE
  126:       ELSE
  127: *
  128: *        Lower triangular storage: examine D from top to bottom.
  129: *
  130:          IP = 1
  131:          DO 20 I = 1, N
  132:             IF( IPIV( I ).GT.0 .AND. AP( IP ).EQ.ZERO )
  133:      $         RETURN
  134:             IP = IP + N - I + 1
  135:    20    CONTINUE
  136:       END IF
  137: *
  138: *     Estimate the 1-norm of the inverse.
  139: *
  140:       KASE = 0
  141:    30 CONTINUE
  142:       CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
  143:       IF( KASE.NE.0 ) THEN
  144: *
  145: *        Multiply by inv(L*D*L') or inv(U*D*U').
  146: *
  147:          CALL ZSPTRS( UPLO, N, 1, AP, IPIV, WORK, N, INFO )
  148:          GO TO 30
  149:       END IF
  150: *
  151: *     Compute the estimate of the reciprocal condition number.
  152: *
  153:       IF( AINVNM.NE.ZERO )
  154:      $   RCOND = ( ONE / AINVNM ) / ANORM
  155: *
  156:       RETURN
  157: *
  158: *     End of ZSPCON
  159: *
  160:       END