Return to zgesc2.f CVS log

Up to [local] / rpl / lapack / lapack

Annotation of rpl/lapack/lapack/zgesc2.f, revision 1.5

1.1       bertrand    1:       SUBROUTINE ZGESC2( N, A, LDA, RHS, IPIV, JPIV, SCALE )
                      2: *
                      3: *  -- LAPACK auxiliary routine (version 3.2) --
                      4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                      5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                      6: *     November 2006
                      7: *
                      8: *     .. Scalar Arguments ..
                      9:       INTEGER            LDA, N
                     10:       DOUBLE PRECISION   SCALE
                     11: *     ..
                     12: *     .. Array Arguments ..
                     13:       INTEGER            IPIV( * ), JPIV( * )
                     14:       COMPLEX*16         A( LDA, * ), RHS( * )
                     15: *     ..
                     16: *
                     17: *  Purpose
                     18: *  =======
                     19: *
                     20: *  ZGESC2 solves a system of linear equations
                     21: *
                     22: *            A * X = scale* RHS
                     23: *
                     24: *  with a general N-by-N matrix A using the LU factorization with
                     25: *  complete pivoting computed by ZGETC2.
                     26: *
                     27: *
                     28: *  Arguments
                     29: *  =========
                     30: *
                     31: *  N       (input) INTEGER
                     32: *          The number of columns of the matrix A.
                     33: *
                     34: *  A       (input) COMPLEX*16 array, dimension (LDA, N)
                     35: *          On entry, the  LU part of the factorization of the n-by-n
                     36: *          matrix A computed by ZGETC2:  A = P * L * U * Q
                     37: *
                     38: *  LDA     (input) INTEGER
                     39: *          The leading dimension of the array A.  LDA >= max(1, N).
                     40: *
                     41: *  RHS     (input/output) COMPLEX*16 array, dimension N.
                     42: *          On entry, the right hand side vector b.
                     43: *          On exit, the solution vector X.
                     44: *
                     45: *  IPIV    (input) INTEGER array, dimension (N).
                     46: *          The pivot indices; for 1 <= i <= N, row i of the
                     47: *          matrix has been interchanged with row IPIV(i).
                     48: *
                     49: *  JPIV    (input) INTEGER array, dimension (N).
                     50: *          The pivot indices; for 1 <= j <= N, column j of the
                     51: *          matrix has been interchanged with column JPIV(j).
                     52: *
                     53: *  SCALE    (output) DOUBLE PRECISION
                     54: *           On exit, SCALE contains the scale factor. SCALE is chosen
                     55: *           0 <= SCALE <= 1 to prevent owerflow in the solution.
                     56: *
                     57: *  Further Details
                     58: *  ===============
                     59: *
                     60: *  Based on contributions by
                     61: *     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
                     62: *     Umea University, S-901 87 Umea, Sweden.
                     63: *
                     64: *  =====================================================================
                     65: *
                     66: *     .. Parameters ..
                     67:       DOUBLE PRECISION   ZERO, ONE, TWO
                     68:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
                     69: *     ..
                     70: *     .. Local Scalars ..
                     71:       INTEGER            I, J
                     72:       DOUBLE PRECISION   BIGNUM, EPS, SMLNUM
                     73:       COMPLEX*16         TEMP
                     74: *     ..
                     75: *     .. External Subroutines ..
                     76:       EXTERNAL           ZLASWP, ZSCAL
                     77: *     ..
                     78: *     .. External Functions ..
                     79:       INTEGER            IZAMAX
                     80:       DOUBLE PRECISION   DLAMCH
                     81:       EXTERNAL           IZAMAX, DLAMCH
                     82: *     ..
                     83: *     .. Intrinsic Functions ..
                     84:       INTRINSIC          ABS, DBLE, DCMPLX
                     85: *     ..
                     86: *     .. Executable Statements ..
                     87: *
                     88: *     Set constant to control overflow
                     89: *
                     90:       EPS = DLAMCH( 'P' )
                     91:       SMLNUM = DLAMCH( 'S' ) / EPS
                     92:       BIGNUM = ONE / SMLNUM
                     93:       CALL DLABAD( SMLNUM, BIGNUM )
                     94: *
                     95: *     Apply permutations IPIV to RHS
                     96: *
                     97:       CALL ZLASWP( 1, RHS, LDA, 1, N-1, IPIV, 1 )
                     98: *
                     99: *     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