001:       SUBROUTINE DGESC2( N, A, LDA, RHS, IPIV, JPIV, SCALE )
002: *
003: *  -- LAPACK auxiliary routine (version 3.2) --
004: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
005: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       INTEGER            LDA, N
010:       DOUBLE PRECISION   SCALE
011: *     ..
012: *     .. Array Arguments ..
013:       INTEGER            IPIV( * ), JPIV( * )
014:       DOUBLE PRECISION   A( LDA, * ), RHS( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  DGESC2 solves a system of linear equations
021: *
022: *            A * X = scale* RHS
023: *
024: *  with a general N-by-N matrix A using the LU factorization with
025: *  complete pivoting computed by DGETC2.
026: *
027: *  Arguments
028: *  =========
029: *
030: *  N       (input) INTEGER
031: *          The order of the matrix A.
032: *
033: *  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
034: *          On entry, the  LU part of the factorization of the n-by-n
035: *          matrix A computed by DGETC2:  A = P * L * U * Q
036: *
037: *  LDA     (input) INTEGER
038: *          The leading dimension of the array A.  LDA >= max(1, N).
039: *
040: *  RHS     (input/output) DOUBLE PRECISION array, dimension (N).
041: *          On entry, the right hand side vector b.
042: *          On exit, the solution vector X.
043: *
044: *  IPIV    (input) INTEGER array, dimension (N).
045: *          The pivot indices; for 1 <= i <= N, row i of the
046: *          matrix has been interchanged with row IPIV(i).
047: *
048: *  JPIV    (input) INTEGER array, dimension (N).
049: *          The pivot indices; for 1 <= j <= N, column j of the
050: *          matrix has been interchanged with column JPIV(j).
051: *
052: *  SCALE    (output) DOUBLE PRECISION
053: *           On exit, SCALE contains the scale factor. SCALE is chosen
054: *           0 <= SCALE <= 1 to prevent owerflow in the solution.
055: *
056: *  Further Details
057: *  ===============
058: *
059: *  Based on contributions by
060: *     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
061: *     Umea University, S-901 87 Umea, Sweden.
062: *
063: *  =====================================================================
064: *
065: *     .. Parameters ..
066:       DOUBLE PRECISION   ONE, TWO
067:       PARAMETER          ( ONE = 1.0D+0, TWO = 2.0D+0 )
068: *     ..
069: *     .. Local Scalars ..
070:       INTEGER            I, J
071:       DOUBLE PRECISION   BIGNUM, EPS, SMLNUM, TEMP
072: *     ..
073: *     .. External Subroutines ..
074:       EXTERNAL           DLASWP, DSCAL
075: *     ..
076: *     .. External Functions ..
077:       INTEGER            IDAMAX
078:       DOUBLE PRECISION   DLAMCH
079:       EXTERNAL           IDAMAX, DLAMCH
080: *     ..
081: *     .. Intrinsic Functions ..
082:       INTRINSIC          ABS
083: *     ..
084: *     .. Executable Statements ..
085: *
086: *      Set constant to control owerflow
087: *
088:       EPS = DLAMCH( 'P' )
089:       SMLNUM = DLAMCH( 'S' ) / EPS
090:       BIGNUM = ONE / SMLNUM
091:       CALL DLABAD( SMLNUM, BIGNUM )
092: *
093: *     Apply permutations IPIV to RHS
094: *
095:       CALL DLASWP( 1, RHS, LDA, 1, N-1, IPIV, 1 )
096: *
097: *     Solve for L part
098: *
099:       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
134: