LAPACK 3.3.0

zgesc2.f

Go to the documentation of this file.
00001       SUBROUTINE ZGESC2( N, A, LDA, RHS, IPIV, JPIV, SCALE )
00002 *
00003 *  -- LAPACK auxiliary routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            LDA, N
00010       DOUBLE PRECISION   SCALE
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            IPIV( * ), JPIV( * )
00014       COMPLEX*16         A( LDA, * ), RHS( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  ZGESC2 solves a system of linear equations
00021 *
00022 *            A * X = scale* RHS
00023 *
00024 *  with a general N-by-N matrix A using the LU factorization with
00025 *  complete pivoting computed by ZGETC2.
00026 *
00027 *
00028 *  Arguments
00029 *  =========
00030 *
00031 *  N       (input) INTEGER
00032 *          The number of columns of the matrix A.
00033 *
00034 *  A       (input) COMPLEX*16 array, dimension (LDA, N)
00035 *          On entry, the  LU part of the factorization of the n-by-n
00036 *          matrix A computed by ZGETC2:  A = P * L * U * Q
00037 *
00038 *  LDA     (input) INTEGER
00039 *          The leading dimension of the array A.  LDA >= max(1, N).
00040 *
00041 *  RHS     (input/output) COMPLEX*16 array, dimension N.
00042 *          On entry, the right hand side vector b.
00043 *          On exit, the solution vector X.
00044 *
00045 *  IPIV    (input) INTEGER array, dimension (N).
00046 *          The pivot indices; for 1 <= i <= N, row i of the
00047 *          matrix has been interchanged with row IPIV(i).
00048 *
00049 *  JPIV    (input) INTEGER array, dimension (N).
00050 *          The pivot indices; for 1 <= j <= N, column j of the
00051 *          matrix has been interchanged with column JPIV(j).
00052 *
00053 *  SCALE    (output) DOUBLE PRECISION
00054 *           On exit, SCALE contains the scale factor. SCALE is chosen
00055 *           0 <= SCALE <= 1 to prevent owerflow in the solution.
00056 *
00057 *  Further Details
00058 *  ===============
00059 *
00060 *  Based on contributions by
00061 *     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
00062 *     Umea University, S-901 87 Umea, Sweden.
00063 *
00064 *  =====================================================================
00065 *
00066 *     .. Parameters ..
00067       DOUBLE PRECISION   ZERO, ONE, TWO
00068       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
00069 *     ..
00070 *     .. Local Scalars ..
00071       INTEGER            I, J
00072       DOUBLE PRECISION   BIGNUM, EPS, SMLNUM
00073       COMPLEX*16         TEMP
00074 *     ..
00075 *     .. External Subroutines ..
00076       EXTERNAL           ZLASWP, ZSCAL
00077 *     ..
00078 *     .. External Functions ..
00079       INTEGER            IZAMAX
00080       DOUBLE PRECISION   DLAMCH
00081       EXTERNAL           IZAMAX, DLAMCH
00082 *     ..
00083 *     .. Intrinsic Functions ..
00084       INTRINSIC          ABS, DBLE, DCMPLX
00085 *     ..
00086 *     .. Executable Statements ..
00087 *
00088 *     Set constant to control overflow
00089 *
00090       EPS = DLAMCH( 'P' )
00091       SMLNUM = DLAMCH( 'S' ) / EPS
00092       BIGNUM = ONE / SMLNUM
00093       CALL DLABAD( SMLNUM, BIGNUM )
00094 *
00095 *     Apply permutations IPIV to RHS
00096 *
00097       CALL ZLASWP( 1, RHS, LDA, 1, N-1, IPIV, 1 )
00098 *
00099 *     Solve for L part
00100 *
00101       DO 20 I = 1, N - 1
00102          DO 10 J = I + 1, N
00103             RHS( J ) = RHS( J ) - A( J, I )*RHS( I )
00104    10    CONTINUE
00105    20 CONTINUE
00106 *
00107 *     Solve for U part
00108 *
00109       SCALE = ONE
00110 *
00111 *     Check for scaling
00112 *
00113       I = IZAMAX( N, RHS, 1 )
00114       IF( TWO*SMLNUM*ABS( RHS( I ) ).GT.ABS( A( N, N ) ) ) THEN
00115          TEMP = DCMPLX( ONE / TWO, ZERO ) / ABS( RHS( I ) )
00116          CALL ZSCAL( N, TEMP, RHS( 1 ), 1 )
00117          SCALE = SCALE*DBLE( TEMP )
00118       END IF
00119       DO 40 I = N, 1, -1
00120          TEMP = DCMPLX( ONE, ZERO ) / A( I, I )
00121          RHS( I ) = RHS( I )*TEMP
00122          DO 30 J = I + 1, N
00123             RHS( I ) = RHS( I ) - RHS( J )*( A( I, J )*TEMP )
00124    30    CONTINUE
00125    40 CONTINUE
00126 *
00127 *     Apply permutations JPIV to the solution (RHS)
00128 *
00129       CALL ZLASWP( 1, RHS, LDA, 1, N-1, JPIV, -1 )
00130       RETURN
00131 *
00132 *     End of ZGESC2
00133 *
00134       END
 All Files Functions