Return to dgesc2.f CVS log

Up to [local] / rpl / lapack / lapack

Annotation of rpl/lapack/lapack/dgesc2.f, revision 1.1.1.1

1.1       bertrand    1:       SUBROUTINE DGESC2( 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:       DOUBLE PRECISION   A( LDA, * ), RHS( * )
                     15: *     ..
                     16: *
                     17: *  Purpose
                     18: *  =======
                     19: *
                     20: *  DGESC2 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 DGETC2.
                     26: *
                     27: *  Arguments
                     28: *  =========
                     29: *
                     30: *  N       (input) INTEGER
                     31: *          The order of the matrix A.
                     32: *
                     33: *  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
                     34: *          On entry, the  LU part of the factorization of the n-by-n
                     35: *          matrix A computed by DGETC2:  A = P * L * U * Q
                     36: *
                     37: *  LDA     (input) INTEGER
                     38: *          The leading dimension of the array A.  LDA >= max(1, N).
                     39: *
                     40: *  RHS     (input/output) DOUBLE PRECISION array, dimension (N).
                     41: *          On entry, the right hand side vector b.
                     42: *          On exit, the solution vector X.
                     43: *
                     44: *  IPIV    (input) INTEGER array, dimension (N).
                     45: *          The pivot indices; for 1 <= i <= N, row i of the
                     46: *          matrix has been interchanged with row IPIV(i).
                     47: *
                     48: *  JPIV    (input) INTEGER array, dimension (N).
                     49: *          The pivot indices; for 1 <= j <= N, column j of the
                     50: *          matrix has been interchanged with column JPIV(j).
                     51: *
                     52: *  SCALE    (output) DOUBLE PRECISION
                     53: *           On exit, SCALE contains the scale factor. SCALE is chosen
                     54: *           0 <= SCALE <= 1 to prevent owerflow in the solution.
                     55: *
                     56: *  Further Details
                     57: *  ===============
                     58: *
                     59: *  Based on contributions by
                     60: *     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
                     61: *     Umea University, S-901 87 Umea, Sweden.
                     62: *
                     63: *  =====================================================================
                     64: *
                     65: *     .. Parameters ..
                     66:       DOUBLE PRECISION   ONE, TWO
                     67:       PARAMETER          ( ONE = 1.0D+0, TWO = 2.0D+0 )
                     68: *     ..
                     69: *     .. Local Scalars ..
                     70:       INTEGER            I, J
                     71:       DOUBLE PRECISION   BIGNUM, EPS, SMLNUM, TEMP
                     72: *     ..
                     73: *     .. External Subroutines ..
                     74:       EXTERNAL           DLASWP, DSCAL
                     75: *     ..
                     76: *     .. External Functions ..
                     77:       INTEGER            IDAMAX
                     78:       DOUBLE PRECISION   DLAMCH
                     79:       EXTERNAL           IDAMAX, DLAMCH
                     80: *     ..
                     81: *     .. Intrinsic Functions ..
                     82:       INTRINSIC          ABS
                     83: *     ..
                     84: *     .. Executable Statements ..
                     85: *
                     86: *      Set constant to control owerflow
                     87: *
                     88:       EPS = DLAMCH( 'P' )
                     89:       SMLNUM = DLAMCH( 'S' ) / EPS
                     90:       BIGNUM = ONE / SMLNUM
                     91:       CALL DLABAD( SMLNUM, BIGNUM )
                     92: *
                     93: *     Apply permutations IPIV to RHS
                     94: *
                     95:       CALL DLASWP( 1, RHS, LDA, 1, N-1, IPIV, 1 )
                     96: *
                     97: *     Solve for L part
                     98: *
                     99:       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