Return to dgecon.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

    1:       SUBROUTINE DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK,
    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 DLACN2 in place of DLACON, 5 Feb 03, SJH.
   10: *
   11: *     .. Scalar Arguments ..
   12:       CHARACTER          NORM
   13:       INTEGER            INFO, LDA, N
   14:       DOUBLE PRECISION   ANORM, RCOND
   15: *     ..
   16: *     .. Array Arguments ..
   17:       INTEGER            IWORK( * )
   18:       DOUBLE PRECISION   A( LDA, * ), WORK( * )
   19: *     ..
   20: *
   21: *  Purpose
   22: *  =======
   23: *
   24: *  DGECON estimates the reciprocal of the condition number of a general
   25: *  real matrix A, in either the 1-norm or the infinity-norm, using
   26: *  the LU factorization computed by DGETRF.
   27: *
   28: *  An estimate is obtained for norm(inv(A)), and the reciprocal of the
   29: *  condition number is computed as
   30: *     RCOND = 1 / ( norm(A) * norm(inv(A)) ).
   31: *
   32: *  Arguments
   33: *  =========
   34: *
   35: *  NORM    (input) CHARACTER*1
   36: *          Specifies whether the 1-norm condition number or the
   37: *          infinity-norm condition number is required:
   38: *          = '1' or 'O':  1-norm;
   39: *          = 'I':         Infinity-norm.
   40: *
   41: *  N       (input) INTEGER
   42: *          The order of the matrix A.  N >= 0.
   43: *
   44: *  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
   45: *          The factors L and U from the factorization A = P*L*U
   46: *          as computed by DGETRF.
   47: *
   48: *  LDA     (input) INTEGER
   49: *          The leading dimension of the array A.  LDA >= max(1,N).
   50: *
   51: *  ANORM   (input) DOUBLE PRECISION
   52: *          If NORM = '1' or 'O', the 1-norm of the original matrix A.
   53: *          If NORM = 'I', the infinity-norm of the original matrix A.
   54: *
   55: *  RCOND   (output) DOUBLE PRECISION
   56: *          The reciprocal of the condition number of the matrix A,
   57: *          computed as RCOND = 1/(norm(A) * norm(inv(A))).
   58: *
   59: *  WORK    (workspace) DOUBLE PRECISION array, dimension (4*N)
   60: *
   61: *  IWORK   (workspace) INTEGER array, dimension (N)
   62: *
   63: *  INFO    (output) INTEGER
   64: *          = 0:  successful exit
   65: *          < 0:  if INFO = -i, the i-th argument had an illegal value
   66: *
   67: *  =====================================================================
   68: *
   69: *     .. Parameters ..
   70:       DOUBLE PRECISION   ONE, ZERO
   71:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
   72: *     ..
   73: *     .. Local Scalars ..
   74:       LOGICAL            ONENRM
   75:       CHARACTER          NORMIN
   76:       INTEGER            IX, KASE, KASE1
   77:       DOUBLE PRECISION   AINVNM, SCALE, SL, SMLNUM, SU
   78: *     ..
   79: *     .. Local Arrays ..
   80:       INTEGER            ISAVE( 3 )
   81: *     ..
   82: *     .. External Functions ..
   83:       LOGICAL            LSAME
   84:       INTEGER            IDAMAX
   85:       DOUBLE PRECISION   DLAMCH
   86:       EXTERNAL           LSAME, IDAMAX, DLAMCH
   87: *     ..
   88: *     .. External Subroutines ..
   89:       EXTERNAL           DLACN2, DLATRS, DRSCL, XERBLA
   90: *     ..
   91: *     .. Intrinsic Functions ..
   92:       INTRINSIC          ABS, MAX
   93: *     ..
   94: *     .. Executable Statements ..
   95: *
   96: *     Test the input parameters.
   97: *
   98:       INFO = 0
   99:       ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
  100:       IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) 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 = -5
  108:       END IF
  109:       IF( INFO.NE.0 ) THEN
  110:          CALL XERBLA( 'DGECON', -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.EQ.ZERO ) THEN
  121:          RETURN
  122:       END IF
  123: *
  124:       SMLNUM = DLAMCH( 'Safe minimum' )
  125: *
  126: *     Estimate the norm of inv(A).
  127: *
  128:       AINVNM = ZERO
  129:       NORMIN = 'N'
  130:       IF( ONENRM ) THEN
  131:          KASE1 = 1
  132:       ELSE
  133:          KASE1 = 2
  134:       END IF
  135:       KASE = 0
  136:    10 CONTINUE
  137:       CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
  138:       IF( KASE.NE.0 ) THEN
  139:          IF( KASE.EQ.KASE1 ) THEN
  140: *
  141: *           Multiply by inv(L).
  142: *
  143:             CALL DLATRS( 'Lower', 'No transpose', 'Unit', NORMIN, N, A,
  144:      $                   LDA, WORK, SL, WORK( 2*N+1 ), INFO )
  145: *
  146: *           Multiply by inv(U).
  147: *
  148:             CALL DLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
  149:      $                   A, LDA, WORK, SU, WORK( 3*N+1 ), INFO )
  150:          ELSE
  151: *
  152: *           Multiply by inv(U').
  153: *
  154:             CALL DLATRS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N, A,
  155:      $                   LDA, WORK, SU, WORK( 3*N+1 ), INFO )
  156: *
  157: *           Multiply by inv(L').
  158: *
  159:             CALL DLATRS( 'Lower', 'Transpose', 'Unit', NORMIN, N, A,
  160:      $                   LDA, WORK, SL, WORK( 2*N+1 ), INFO )
  161:          END IF
  162: *
  163: *        Divide X by 1/(SL*SU) if doing so will not cause overflow.
  164: *
  165:          SCALE = SL*SU
  166:          NORMIN = 'Y'
  167:          IF( SCALE.NE.ONE ) THEN
  168:             IX = IDAMAX( N, WORK, 1 )
  169:             IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
  170:      $         GO TO 20
  171:             CALL DRSCL( N, SCALE, WORK, 1 )
  172:          END IF
  173:          GO TO 10
  174:       END IF
  175: *
  176: *     Compute the estimate of the reciprocal condition number.
  177: *
  178:       IF( AINVNM.NE.ZERO )
  179:      $   RCOND = ( ONE / AINVNM ) / ANORM
  180: *
  181:    20 CONTINUE
  182:       RETURN
  183: *
  184: *     End of DGECON
  185: *
  186:       END