File:  [local] / rpl / lapack / lapack / dgesc2.f
Revision 1.21: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:38:49 2023 UTC (8 months, 3 weeks ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_35, rpl-4_1_34, HEAD
Première mise à jour de lapack et blas.

    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 overflow 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: *> \ingroup doubleGEauxiliary
  105: *
  106: *> \par Contributors:
  107: *  ==================
  108: *>
  109: *>     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
  110: *>     Umea University, S-901 87 Umea, Sweden.
  111: *
  112: *  =====================================================================
  113:       SUBROUTINE DGESC2( N, A, LDA, RHS, IPIV, JPIV, SCALE )
  114: *
  115: *  -- LAPACK auxiliary routine --
  116: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  117: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  118: *
  119: *     .. Scalar Arguments ..
  120:       INTEGER            LDA, N
  121:       DOUBLE PRECISION   SCALE
  122: *     ..
  123: *     .. Array Arguments ..
  124:       INTEGER            IPIV( * ), JPIV( * )
  125:       DOUBLE PRECISION   A( LDA, * ), RHS( * )
  126: *     ..
  127: *
  128: *  =====================================================================
  129: *
  130: *     .. Parameters ..
  131:       DOUBLE PRECISION   ONE, TWO
  132:       PARAMETER          ( ONE = 1.0D+0, TWO = 2.0D+0 )
  133: *     ..
  134: *     .. Local Scalars ..
  135:       INTEGER            I, J
  136:       DOUBLE PRECISION   BIGNUM, EPS, SMLNUM, TEMP
  137: *     ..
  138: *     .. External Subroutines ..
  139:       EXTERNAL           DLASWP, DSCAL, DLABAD
  140: *     ..
  141: *     .. External Functions ..
  142:       INTEGER            IDAMAX
  143:       DOUBLE PRECISION   DLAMCH
  144:       EXTERNAL           IDAMAX, DLAMCH
  145: *     ..
  146: *     .. Intrinsic Functions ..
  147:       INTRINSIC          ABS
  148: *     ..
  149: *     .. Executable Statements ..
  150: *
  151: *      Set constant to control overflow
  152: *
  153:       EPS = DLAMCH( 'P' )
  154:       SMLNUM = DLAMCH( 'S' ) / EPS
  155:       BIGNUM = ONE / SMLNUM
  156:       CALL DLABAD( SMLNUM, BIGNUM )
  157: *
  158: *     Apply permutations IPIV to RHS
  159: *
  160:       CALL DLASWP( 1, RHS, LDA, 1, N-1, IPIV, 1 )
  161: *
  162: *     Solve for L part
  163: *
  164:       DO 20 I = 1, N - 1
  165:          DO 10 J = I + 1, N
  166:             RHS( J ) = RHS( J ) - A( J, I )*RHS( I )
  167:    10    CONTINUE
  168:    20 CONTINUE
  169: *
  170: *     Solve for U part
  171: *
  172:       SCALE = ONE
  173: *
  174: *     Check for scaling
  175: *
  176:       I = IDAMAX( N, RHS, 1 )
  177:       IF( TWO*SMLNUM*ABS( RHS( I ) ).GT.ABS( A( N, N ) ) ) THEN
  178:          TEMP = ( ONE / TWO ) / ABS( RHS( I ) )
  179:          CALL DSCAL( N, TEMP, RHS( 1 ), 1 )
  180:          SCALE = SCALE*TEMP
  181:       END IF
  182: *
  183:       DO 40 I = N, 1, -1
  184:          TEMP = ONE / A( I, I )
  185:          RHS( I ) = RHS( I )*TEMP
  186:          DO 30 J = I + 1, N
  187:             RHS( I ) = RHS( I ) - RHS( J )*( A( I, J )*TEMP )
  188:    30    CONTINUE
  189:    40 CONTINUE
  190: *
  191: *     Apply permutations JPIV to the solution (RHS)
  192: *
  193:       CALL DLASWP( 1, RHS, LDA, 1, N-1, JPIV, -1 )
  194:       RETURN
  195: *
  196: *     End of DGESC2
  197: *
  198:       END

CVSweb interface <joel.bertrand@systella.fr>