Return to dgesc2.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de Lapack.

    1: *> \brief \b DGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed by sgetc2.
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download DGESC2 + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesc2.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesc2.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesc2.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DGESC2( N, A, LDA, RHS, IPIV, JPIV, SCALE )
   22: *
   23: *       .. Scalar Arguments ..
   24: *       INTEGER            LDA, N
   25: *       DOUBLE PRECISION   SCALE
   26: *       ..
   27: *       .. Array Arguments ..
   28: *       INTEGER            IPIV( * ), JPIV( * )
   29: *       DOUBLE PRECISION   A( LDA, * ), RHS( * )
   30: *       ..
   31: *
   32: *
   33: *> \par Purpose:
   34: *  =============
   35: *>
   36: *> \verbatim
   37: *>
   38: *> DGESC2 solves a system of linear equations
   39: *>
   40: *>           A * X = scale* RHS
   41: *>
   42: *> with a general N-by-N matrix A using the LU factorization with
   43: *> complete pivoting computed by DGETC2.
   44: *> \endverbatim
   45: *
   46: *  Arguments:
   47: *  ==========
   48: *
   49: *> \param[in] N
   50: *> \verbatim
   51: *>          N is INTEGER
   52: *>          The order of the matrix A.
   53: *> \endverbatim
   54: *>
   55: *> \param[in] A
   56: *> \verbatim
   57: *>          A is DOUBLE PRECISION array, dimension (LDA,N)
   58: *>          On entry, the  LU part of the factorization of the n-by-n
   59: *>          matrix A computed by DGETC2:  A = P * L * U * Q
   60: *> \endverbatim
   61: *>
   62: *> \param[in] LDA
   63: *> \verbatim
   64: *>          LDA is INTEGER
   65: *>          The leading dimension of the array A.  LDA >= max(1, N).
   66: *> \endverbatim
   67: *>
   68: *> \param[in,out] RHS
   69: *> \verbatim
   70: *>          RHS is DOUBLE PRECISION array, dimension (N).
   71: *>          On entry, the right hand side vector b.
   72: *>          On exit, the solution vector X.
   73: *> \endverbatim
   74: *>
   75: *> \param[in] IPIV
   76: *> \verbatim
   77: *>          IPIV is INTEGER array, dimension (N).
   78: *>          The pivot indices; for 1 <= i <= N, row i of the
   79: *>          matrix has been interchanged with row IPIV(i).
   80: *> \endverbatim
   81: *>
   82: *> \param[in] JPIV
   83: *> \verbatim
   84: *>          JPIV is INTEGER array, dimension (N).
   85: *>          The pivot indices; for 1 <= j <= N, column j of the
   86: *>          matrix has been interchanged with column JPIV(j).
   87: *> \endverbatim
   88: *>
   89: *> \param[out] SCALE
   90: *> \verbatim
   91: *>          SCALE is DOUBLE PRECISION
   92: *>          On exit, SCALE contains the scale factor. SCALE is chosen
   93: *>          0 <= SCALE <= 1 to prevent owerflow in the solution.
   94: *> \endverbatim
   95: *
   96: *  Authors:
   97: *  ========
   98: *
   99: *> \author Univ. of Tennessee
  100: *> \author Univ. of California Berkeley
  101: *> \author Univ. of Colorado Denver
  102: *> \author NAG Ltd.
  103: *
  104: *> \date November 2017
  105: *
  106: *> \ingroup doubleGEauxiliary
  107: *
  108: *> \par Contributors:
  109: *  ==================
  110: *>
  111: *>     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
  112: *>     Umea University, S-901 87 Umea, Sweden.
  113: *
  114: *  =====================================================================
  115:       SUBROUTINE DGESC2( N, A, LDA, RHS, IPIV, JPIV, SCALE )
  116: *
  117: *  -- LAPACK auxiliary routine (version 3.8.0) --
  118: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  119: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  120: *     November 2017
  121: *
  122: *     .. Scalar Arguments ..
  123:       INTEGER            LDA, N
  124:       DOUBLE PRECISION   SCALE
  125: *     ..
  126: *     .. Array Arguments ..
  127:       INTEGER            IPIV( * ), JPIV( * )
  128:       DOUBLE PRECISION   A( LDA, * ), RHS( * )
  129: *     ..
  130: *
  131: *  =====================================================================
  132: *
  133: *     .. Parameters ..
  134:       DOUBLE PRECISION   ONE, TWO
  135:       PARAMETER          ( ONE = 1.0D+0, TWO = 2.0D+0 )
  136: *     ..
  137: *     .. Local Scalars ..
  138:       INTEGER            I, J
  139:       DOUBLE PRECISION   BIGNUM, EPS, SMLNUM, TEMP
  140: *     ..
  141: *     .. External Subroutines ..
  142:       EXTERNAL           DLASWP, DSCAL, DLABAD
  143: *     ..
  144: *     .. External Functions ..
  145:       INTEGER            IDAMAX
  146:       DOUBLE PRECISION   DLAMCH
  147:       EXTERNAL           IDAMAX, DLAMCH
  148: *     ..
  149: *     .. Intrinsic Functions ..
  150:       INTRINSIC          ABS
  151: *     ..
  152: *     .. Executable Statements ..
  153: *
  154: *      Set constant to control owerflow
  155: *
  156:       EPS = DLAMCH( 'P' )
  157:       SMLNUM = DLAMCH( 'S' ) / EPS
  158:       BIGNUM = ONE / SMLNUM
  159:       CALL DLABAD( SMLNUM, BIGNUM )
  160: *
  161: *     Apply permutations IPIV to RHS
  162: *
  163:       CALL DLASWP( 1, RHS, LDA, 1, N-1, IPIV, 1 )
  164: *
  165: *     Solve for L part
  166: *
  167:       DO 20 I = 1, N - 1
  168:          DO 10 J = I + 1, N
  169:             RHS( J ) = RHS( J ) - A( J, I )*RHS( I )
  170:    10    CONTINUE
  171:    20 CONTINUE
  172: *
  173: *     Solve for U part
  174: *
  175:       SCALE = ONE
  176: *
  177: *     Check for scaling
  178: *
  179:       I = IDAMAX( N, RHS, 1 )
  180:       IF( TWO*SMLNUM*ABS( RHS( I ) ).GT.ABS( A( N, N ) ) ) THEN
  181:          TEMP = ( ONE / TWO ) / ABS( RHS( I ) )
  182:          CALL DSCAL( N, TEMP, RHS( 1 ), 1 )
  183:          SCALE = SCALE*TEMP
  184:       END IF
  185: *
  186:       DO 40 I = N, 1, -1
  187:          TEMP = ONE / A( I, I )
  188:          RHS( I ) = RHS( I )*TEMP
  189:          DO 30 J = I + 1, N
  190:             RHS( I ) = RHS( I ) - RHS( J )*( A( I, J )*TEMP )
  191:    30    CONTINUE
  192:    40 CONTINUE
  193: *
  194: *     Apply permutations JPIV to the solution (RHS)
  195: *
  196:       CALL DLASWP( 1, RHS, LDA, 1, N-1, JPIV, -1 )
  197:       RETURN
  198: *
  199: *     End of DGESC2
  200: *
  201:       END