001:       SUBROUTINE ZGESC2( 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:       COMPLEX*16         A( LDA, * ), RHS( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  ZGESC2 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 ZGETC2.
026: *
027: *
028: *  Arguments
029: *  =========
030: *
031: *  N       (input) INTEGER
032: *          The number of columns of the matrix A.
033: *
034: *  A       (input) COMPLEX*16 array, dimension (LDA, N)
035: *          On entry, the  LU part of the factorization of the n-by-n
036: *          matrix A computed by ZGETC2:  A = P * L * U * Q
037: *
038: *  LDA     (input) INTEGER
039: *          The leading dimension of the array A.  LDA >= max(1, N).
040: *
041: *  RHS     (input/output) COMPLEX*16 array, dimension N.
042: *          On entry, the right hand side vector b.
043: *          On exit, the solution vector X.
044: *
045: *  IPIV    (input) INTEGER array, dimension (N).
046: *          The pivot indices; for 1 <= i <= N, row i of the
047: *          matrix has been interchanged with row IPIV(i).
048: *
049: *  JPIV    (input) INTEGER array, dimension (N).
050: *          The pivot indices; for 1 <= j <= N, column j of the
051: *          matrix has been interchanged with column JPIV(j).
052: *
053: *  SCALE    (output) DOUBLE PRECISION
054: *           On exit, SCALE contains the scale factor. SCALE is chosen
055: *           0 <= SCALE <= 1 to prevent owerflow in the solution.
056: *
057: *  Further Details
058: *  ===============
059: *
060: *  Based on contributions by
061: *     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
062: *     Umea University, S-901 87 Umea, Sweden.
063: *
064: *  =====================================================================
065: *
066: *     .. Parameters ..
067:       DOUBLE PRECISION   ZERO, ONE, TWO
068:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
069: *     ..
070: *     .. Local Scalars ..
071:       INTEGER            I, J
072:       DOUBLE PRECISION   BIGNUM, EPS, SMLNUM
073:       COMPLEX*16         TEMP
074: *     ..
075: *     .. External Subroutines ..
076:       EXTERNAL           ZLASWP, ZSCAL
077: *     ..
078: *     .. External Functions ..
079:       INTEGER            IZAMAX
080:       DOUBLE PRECISION   DLAMCH
081:       EXTERNAL           IZAMAX, DLAMCH
082: *     ..
083: *     .. Intrinsic Functions ..
084:       INTRINSIC          ABS, DBLE, DCMPLX
085: *     ..
086: *     .. Executable Statements ..
087: *
088: *     Set constant to control overflow
089: *
090:       EPS = DLAMCH( 'P' )
091:       SMLNUM = DLAMCH( 'S' ) / EPS
092:       BIGNUM = ONE / SMLNUM
093:       CALL DLABAD( SMLNUM, BIGNUM )
094: *
095: *     Apply permutations IPIV to RHS
096: *
097:       CALL ZLASWP( 1, RHS, LDA, 1, N-1, IPIV, 1 )
098: *
099: *     Solve for L part
100: *
101:       DO 20 I = 1, N - 1
102:          DO 10 J = I + 1, N
103:             RHS( J ) = RHS( J ) - A( J, I )*RHS( I )
104:    10    CONTINUE
105:    20 CONTINUE
106: *
107: *     Solve for U part
108: *
109:       SCALE = ONE
110: *
111: *     Check for scaling
112: *
113:       I = IZAMAX( N, RHS, 1 )
114:       IF( TWO*SMLNUM*ABS( RHS( I ) ).GT.ABS( A( N, N ) ) ) THEN
115:          TEMP = DCMPLX( ONE / TWO, ZERO ) / ABS( RHS( I ) )
116:          CALL ZSCAL( N, TEMP, RHS( 1 ), 1 )
117:          SCALE = SCALE*DBLE( TEMP )
118:       END IF
119:       DO 40 I = N, 1, -1
120:          TEMP = DCMPLX( ONE, ZERO ) / A( I, I )
121:          RHS( I ) = RHS( I )*TEMP
122:          DO 30 J = I + 1, N
123:             RHS( I ) = RHS( I ) - RHS( J )*( A( I, J )*TEMP )
124:    30    CONTINUE
125:    40 CONTINUE
126: *
127: *     Apply permutations JPIV to the solution (RHS)
128: *
129:       CALL ZLASWP( 1, RHS, LDA, 1, N-1, JPIV, -1 )
130:       RETURN
131: *
132: *     End of ZGESC2
133: *
134:       END
135: