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