Return to dgesc2.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour globale de Lapack 3.2.2.

    1:       SUBROUTINE DGESC2( N, A, LDA, RHS, IPIV, JPIV, SCALE )
    2: *
    3: *  -- LAPACK auxiliary routine (version 3.2.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: *     June 2010
    7: *
    8: *     .. Scalar Arguments ..
    9:       INTEGER            LDA, N
   10:       DOUBLE PRECISION   SCALE
   11: *     ..
   12: *     .. Array Arguments ..
   13:       INTEGER            IPIV( * ), JPIV( * )
   14:       DOUBLE PRECISION   A( LDA, * ), RHS( * )
   15: *     ..
   16: *
   17: *  Purpose
   18: *  =======
   19: *
   20: *  DGESC2 solves a system of linear equations
   21: *
   22: *            A * X = scale* RHS
   23: *
   24: *  with a general N-by-N matrix A using the LU factorization with
   25: *  complete pivoting computed by DGETC2.
   26: *
   27: *  Arguments
   28: *  =========
   29: *
   30: *  N       (input) INTEGER
   31: *          The order of the matrix A.
   32: *
   33: *  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
   34: *          On entry, the  LU part of the factorization of the n-by-n
   35: *          matrix A computed by DGETC2:  A = P * L * U * Q
   36: *
   37: *  LDA     (input) INTEGER
   38: *          The leading dimension of the array A.  LDA >= max(1, N).
   39: *
   40: *  RHS     (input/output) DOUBLE PRECISION array, dimension (N).
   41: *          On entry, the right hand side vector b.
   42: *          On exit, the solution vector X.
   43: *
   44: *  IPIV    (input) INTEGER array, dimension (N).
   45: *          The pivot indices; for 1 <= i <= N, row i of the
   46: *          matrix has been interchanged with row IPIV(i).
   47: *
   48: *  JPIV    (input) INTEGER array, dimension (N).
   49: *          The pivot indices; for 1 <= j <= N, column j of the
   50: *          matrix has been interchanged with column JPIV(j).
   51: *
   52: *  SCALE   (output) DOUBLE PRECISION
   53: *          On exit, SCALE contains the scale factor. SCALE is chosen
   54: *          0 <= SCALE <= 1 to prevent owerflow in the solution.
   55: *
   56: *  Further Details
   57: *  ===============
   58: *
   59: *  Based on contributions by
   60: *     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
   61: *     Umea University, S-901 87 Umea, Sweden.
   62: *
   63: *  =====================================================================
   64: *
   65: *     .. Parameters ..
   66:       DOUBLE PRECISION   ONE, TWO
   67:       PARAMETER          ( ONE = 1.0D+0, TWO = 2.0D+0 )
   68: *     ..
   69: *     .. Local Scalars ..
   70:       INTEGER            I, J
   71:       DOUBLE PRECISION   BIGNUM, EPS, SMLNUM, TEMP
   72: *     ..
   73: *     .. External Subroutines ..
   74:       EXTERNAL           DLASWP, DSCAL
   75: *     ..
   76: *     .. External Functions ..
   77:       INTEGER            IDAMAX
   78:       DOUBLE PRECISION   DLAMCH
   79:       EXTERNAL           IDAMAX, DLAMCH
   80: *     ..
   81: *     .. Intrinsic Functions ..
   82:       INTRINSIC          ABS
   83: *     ..
   84: *     .. Executable Statements ..
   85: *
   86: *      Set constant to control owerflow
   87: *
   88:       EPS = DLAMCH( 'P' )
   89:       SMLNUM = DLAMCH( 'S' ) / EPS
   90:       BIGNUM = ONE / SMLNUM
   91:       CALL DLABAD( SMLNUM, BIGNUM )
   92: *
   93: *     Apply permutations IPIV to RHS
   94: *
   95:       CALL DLASWP( 1, RHS, LDA, 1, N-1, IPIV, 1 )
   96: *
   97: *     Solve for L part
   98: *
   99:       DO 20 I = 1, N - 1
  100:          DO 10 J = I + 1, N
  101:             RHS( J ) = RHS( J ) - A( J, I )*RHS( I )
  102:    10    CONTINUE
  103:    20 CONTINUE
  104: *
  105: *     Solve for U part
  106: *
  107:       SCALE = ONE
  108: *
  109: *     Check for scaling
  110: *
  111:       I = IDAMAX( N, RHS, 1 )
  112:       IF( TWO*SMLNUM*ABS( RHS( I ) ).GT.ABS( A( N, N ) ) ) THEN
  113:          TEMP = ( ONE / TWO ) / ABS( RHS( I ) )
  114:          CALL DSCAL( N, TEMP, RHS( 1 ), 1 )
  115:          SCALE = SCALE*TEMP
  116:       END IF
  117: *
  118:       DO 40 I = N, 1, -1
  119:          TEMP = ONE / A( I, I )
  120:          RHS( I ) = RHS( I )*TEMP
  121:          DO 30 J = I + 1, N
  122:             RHS( I ) = RHS( I ) - RHS( J )*( A( I, J )*TEMP )
  123:    30    CONTINUE
  124:    40 CONTINUE
  125: *
  126: *     Apply permutations JPIV to the solution (RHS)
  127: *
  128:       CALL DLASWP( 1, RHS, LDA, 1, N-1, JPIV, -1 )
  129:       RETURN
  130: *
  131: *     End of DGESC2
  132: *
  133:       END