diff options
Diffstat (limited to '2.3-1/src/fortran/lapack/zgesc2.f')
| -rw-r--r-- | 2.3-1/src/fortran/lapack/zgesc2.f | 133 |
1 files changed, 133 insertions, 0 deletions
diff --git a/2.3-1/src/fortran/lapack/zgesc2.f b/2.3-1/src/fortran/lapack/zgesc2.f new file mode 100644 index 00000000..d4d51337 --- /dev/null +++ b/2.3-1/src/fortran/lapack/zgesc2.f @@ -0,0 +1,133 @@ + SUBROUTINE ZGESC2( N, A, LDA, RHS, IPIV, JPIV, SCALE ) +* +* -- LAPACK auxiliary routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + INTEGER LDA, N + DOUBLE PRECISION SCALE +* .. +* .. Array Arguments .. + INTEGER IPIV( * ), JPIV( * ) + COMPLEX*16 A( LDA, * ), RHS( * ) +* .. +* +* Purpose +* ======= +* +* ZGESC2 solves a system of linear equations +* +* A * X = scale* RHS +* +* with a general N-by-N matrix A using the LU factorization with +* complete pivoting computed by ZGETC2. +* +* +* Arguments +* ========= +* +* N (input) INTEGER +* The number of columns of the matrix A. +* +* A (input) COMPLEX*16 array, dimension (LDA, N) +* On entry, the LU part of the factorization of the n-by-n +* matrix A computed by ZGETC2: A = P * L * U * Q +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1, N). +* +* RHS (input/output) COMPLEX*16 array, dimension N. +* On entry, the right hand side vector b. +* On exit, the solution vector X. +* +* IPIV (input) INTEGER array, dimension (N). +* The pivot indices; for 1 <= i <= N, row i of the +* matrix has been interchanged with row IPIV(i). +* +* JPIV (input) INTEGER array, dimension (N). +* The pivot indices; for 1 <= j <= N, column j of the +* matrix has been interchanged with column JPIV(j). +* +* SCALE (output) DOUBLE PRECISION +* On exit, SCALE contains the scale factor. SCALE is chosen +* 0 <= SCALE <= 1 to prevent owerflow in the solution. +* +* Further Details +* =============== +* +* Based on contributions by +* Bo Kagstrom and Peter Poromaa, Department of Computing Science, +* Umea University, S-901 87 Umea, Sweden. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO, ONE, TWO + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 ) +* .. +* .. Local Scalars .. + INTEGER I, J + DOUBLE PRECISION BIGNUM, EPS, SMLNUM + COMPLEX*16 TEMP +* .. +* .. External Subroutines .. + EXTERNAL ZLASWP, ZSCAL +* .. +* .. External Functions .. + INTEGER IZAMAX + DOUBLE PRECISION DLAMCH + EXTERNAL IZAMAX, DLAMCH +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, DBLE, DCMPLX +* .. +* .. Executable Statements .. +* +* Set constant to control overflow +* + EPS = DLAMCH( 'P' ) + SMLNUM = DLAMCH( 'S' ) / EPS + BIGNUM = ONE / SMLNUM + CALL DLABAD( SMLNUM, BIGNUM ) +* +* Apply permutations IPIV to RHS +* + CALL ZLASWP( 1, RHS, LDA, 1, N-1, IPIV, 1 ) +* +* Solve for L part +* + DO 20 I = 1, N - 1 + DO 10 J = I + 1, N + RHS( J ) = RHS( J ) - A( J, I )*RHS( I ) + 10 CONTINUE + 20 CONTINUE +* +* Solve for U part +* + SCALE = ONE +* +* Check for scaling +* + I = IZAMAX( N, RHS, 1 ) + IF( TWO*SMLNUM*ABS( RHS( I ) ).GT.ABS( A( N, N ) ) ) THEN + TEMP = DCMPLX( ONE / TWO, ZERO ) / ABS( RHS( I ) ) + CALL ZSCAL( N, TEMP, RHS( 1 ), 1 ) + SCALE = SCALE*DBLE( TEMP ) + END IF + DO 40 I = N, 1, -1 + TEMP = DCMPLX( ONE, ZERO ) / A( I, I ) + RHS( I ) = RHS( I )*TEMP + DO 30 J = I + 1, N + RHS( I ) = RHS( I ) - RHS( J )*( A( I, J )*TEMP ) + 30 CONTINUE + 40 CONTINUE +* +* Apply permutations JPIV to the solution (RHS) +* + CALL ZLASWP( 1, RHS, LDA, 1, N-1, JPIV, -1 ) + RETURN +* +* End of ZGESC2 +* + END |