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