Return to zgesc2.f CVS log

Up to [local] / rpl / lapack / lapack

En route pour la 4.0.15 !

    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