diff options
Diffstat (limited to 'SRC/dgesc2.f')
| -rw-r--r-- | SRC/dgesc2.f | 132 |
1 files changed, 132 insertions, 0 deletions
diff --git a/SRC/dgesc2.f b/SRC/dgesc2.f new file mode 100644 index 00000000..1b0331f5 --- /dev/null +++ b/SRC/dgesc2.f @@ -0,0 +1,132 @@ + SUBROUTINE DGESC2( 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( * ) + DOUBLE PRECISION A( LDA, * ), RHS( * ) +* .. +* +* Purpose +* ======= +* +* DGESC2 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 DGETC2. +* +* Arguments +* ========= +* +* N (input) INTEGER +* The order of the matrix A. +* +* A (input) DOUBLE PRECISION array, dimension (LDA,N) +* On entry, the LU part of the factorization of the n-by-n +* matrix A computed by DGETC2: A = P * L * U * Q +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1, N). +* +* RHS (input/output) DOUBLE PRECISION 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 ONE, TWO + PARAMETER ( ONE = 1.0D+0, TWO = 2.0D+0 ) +* .. +* .. Local Scalars .. + INTEGER I, J + DOUBLE PRECISION BIGNUM, EPS, SMLNUM, TEMP +* .. +* .. External Subroutines .. + EXTERNAL DLASWP, DSCAL +* .. +* .. External Functions .. + INTEGER IDAMAX + DOUBLE PRECISION DLAMCH + EXTERNAL IDAMAX, DLAMCH +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS +* .. +* .. Executable Statements .. +* +* Set constant to control owerflow +* + EPS = DLAMCH( 'P' ) + SMLNUM = DLAMCH( 'S' ) / EPS + BIGNUM = ONE / SMLNUM + CALL DLABAD( SMLNUM, BIGNUM ) +* +* Apply permutations IPIV to RHS +* + CALL DLASWP( 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 = IDAMAX( N, RHS, 1 ) + IF( TWO*SMLNUM*ABS( RHS( I ) ).GT.ABS( A( N, N ) ) ) THEN + TEMP = ( ONE / TWO ) / ABS( RHS( I ) ) + CALL DSCAL( N, TEMP, RHS( 1 ), 1 ) + SCALE = SCALE*TEMP + END IF +* + DO 40 I = N, 1, -1 + TEMP = ONE / 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 DLASWP( 1, RHS, LDA, 1, N-1, JPIV, -1 ) + RETURN +* +* End of DGESC2 +* + END |