Return to dsycon.f CVS log

Up to [local] / rpl / lapack / lapack

Patches pour OS/2

    1:       SUBROUTINE DSYCON( UPLO, N, A, LDA, IPIV, ANORM, RCOND, WORK,
    2:      $                   IWORK, INFO )
    3: *
    4: *  -- LAPACK routine (version 3.2.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: *     June 2010
    8: *
    9: *     Modified to call DLACN2 in place of DLACON, 5 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:       INTEGER            IPIV( * ), IWORK( * )
   18:       DOUBLE PRECISION   A( LDA, * ), WORK( * )
   19: *     ..
   20: *
   21: *  Purpose
   22: *  =======
   23: *
   24: *  DSYCON estimates the reciprocal of the condition number (in the
   25: *  1-norm) of a real symmetric matrix A using the factorization
   26: *  A = U*D*U**T or A = L*D*L**T computed by DSYTRF.
   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: *          Specifies whether the details of the factorization are stored
   36: *          as an upper or lower triangular matrix.
   37: *          = 'U':  Upper triangular, form is A = U*D*U**T;
   38: *          = 'L':  Lower triangular, form is A = L*D*L**T.
   39: *
   40: *  N       (input) INTEGER
   41: *          The order of the matrix A.  N >= 0.
   42: *
   43: *  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
   44: *          The block diagonal matrix D and the multipliers used to
   45: *          obtain the factor U or L as computed by DSYTRF.
   46: *
   47: *  LDA     (input) INTEGER
   48: *          The leading dimension of the array A.  LDA >= max(1,N).
   49: *
   50: *  IPIV    (input) INTEGER array, dimension (N)
   51: *          Details of the interchanges and the block structure of D
   52: *          as determined by DSYTRF.
   53: *
   54: *  ANORM   (input) DOUBLE PRECISION
   55: *          The 1-norm of the original matrix A.
   56: *
   57: *  RCOND   (output) DOUBLE PRECISION
   58: *          The reciprocal of the condition number of the matrix A,
   59: *          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
   60: *          estimate of the 1-norm of inv(A) computed in this routine.
   61: *
   62: *  WORK    (workspace) DOUBLE PRECISION array, dimension (2*N)
   63: *
   64: *  IWORK   (workspace) INTEGER array, dimension (N)
   65: *
   66: *  INFO    (output) INTEGER
   67: *          = 0:  successful exit
   68: *          < 0:  if INFO = -i, the i-th argument had an illegal value
   69: *
   70: *  =====================================================================
   71: *
   72: *     .. Parameters ..
   73:       DOUBLE PRECISION   ONE, ZERO
   74:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
   75: *     ..
   76: *     .. Local Scalars ..
   77:       LOGICAL            UPPER
   78:       INTEGER            I, KASE
   79:       DOUBLE PRECISION   AINVNM
   80: *     ..
   81: *     .. Local Arrays ..
   82:       INTEGER            ISAVE( 3 )
   83: *     ..
   84: *     .. External Functions ..
   85:       LOGICAL            LSAME
   86:       EXTERNAL           LSAME
   87: *     ..
   88: *     .. External Subroutines ..
   89:       EXTERNAL           DLACN2, DSYTRS, XERBLA
   90: *     ..
   91: *     .. Intrinsic Functions ..
   92:       INTRINSIC          MAX
   93: *     ..
   94: *     .. Executable Statements ..
   95: *
   96: *     Test the input parameters.
   97: *
   98:       INFO = 0
   99:       UPPER = LSAME( UPLO, 'U' )
  100:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  101:          INFO = -1
  102:       ELSE IF( N.LT.0 ) THEN
  103:          INFO = -2
  104:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  105:          INFO = -4
  106:       ELSE IF( ANORM.LT.ZERO ) THEN
  107:          INFO = -6
  108:       END IF
  109:       IF( INFO.NE.0 ) THEN
  110:          CALL XERBLA( 'DSYCON', -INFO )
  111:          RETURN
  112:       END IF
  113: *
  114: *     Quick return if possible
  115: *
  116:       RCOND = ZERO
  117:       IF( N.EQ.0 ) THEN
  118:          RCOND = ONE
  119:          RETURN
  120:       ELSE IF( ANORM.LE.ZERO ) THEN
  121:          RETURN
  122:       END IF
  123: *
  124: *     Check that the diagonal matrix D is nonsingular.
  125: *
  126:       IF( UPPER ) THEN
  127: *
  128: *        Upper triangular storage: examine D from bottom to top
  129: *
  130:          DO 10 I = N, 1, -1
  131:             IF( IPIV( I ).GT.0 .AND. A( I, I ).EQ.ZERO )
  132:      $         RETURN
  133:    10    CONTINUE
  134:       ELSE
  135: *
  136: *        Lower triangular storage: examine D from top to bottom.
  137: *
  138:          DO 20 I = 1, N
  139:             IF( IPIV( I ).GT.0 .AND. A( I, I ).EQ.ZERO )
  140:      $         RETURN
  141:    20    CONTINUE
  142:       END IF
  143: *
  144: *     Estimate the 1-norm of the inverse.
  145: *
  146:       KASE = 0
  147:    30 CONTINUE
  148:       CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
  149:       IF( KASE.NE.0 ) THEN
  150: *
  151: *        Multiply by inv(L*D*L') or inv(U*D*U').
  152: *
  153:          CALL DSYTRS( UPLO, N, 1, A, LDA, IPIV, WORK, N, INFO )
  154:          GO TO 30
  155:       END IF
  156: *
  157: *     Compute the estimate of the reciprocal condition number.
  158: *
  159:       IF( AINVNM.NE.ZERO )
  160:      $   RCOND = ( ONE / AINVNM ) / ANORM
  161: *
  162:       RETURN
  163: *
  164: *     End of DSYCON
  165: *
  166:       END