Annotation of rpl/lapack/lapack/dgesc2.f, revision 1.3
1.1 bertrand 1: SUBROUTINE DGESC2( 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: DOUBLE PRECISION A( LDA, * ), RHS( * )
15: * ..
16: *
17: * Purpose
18: * =======
19: *
20: * DGESC2 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 DGETC2.
26: *
27: * Arguments
28: * =========
29: *
30: * N (input) INTEGER
31: * The order of the matrix A.
32: *
33: * A (input) DOUBLE PRECISION array, dimension (LDA,N)
34: * On entry, the LU part of the factorization of the n-by-n
35: * matrix A computed by DGETC2: A = P * L * U * Q
36: *
37: * LDA (input) INTEGER
38: * The leading dimension of the array A. LDA >= max(1, N).
39: *
40: * RHS (input/output) DOUBLE PRECISION array, dimension (N).
41: * On entry, the right hand side vector b.
42: * On exit, the solution vector X.
43: *
44: * IPIV (input) INTEGER array, dimension (N).
45: * The pivot indices; for 1 <= i <= N, row i of the
46: * matrix has been interchanged with row IPIV(i).
47: *
48: * JPIV (input) INTEGER array, dimension (N).
49: * The pivot indices; for 1 <= j <= N, column j of the
50: * matrix has been interchanged with column JPIV(j).
51: *
52: * SCALE (output) DOUBLE PRECISION
53: * On exit, SCALE contains the scale factor. SCALE is chosen
54: * 0 <= SCALE <= 1 to prevent owerflow in the solution.
55: *
56: * Further Details
57: * ===============
58: *
59: * Based on contributions by
60: * Bo Kagstrom and Peter Poromaa, Department of Computing Science,
61: * Umea University, S-901 87 Umea, Sweden.
62: *
63: * =====================================================================
64: *
65: * .. Parameters ..
66: DOUBLE PRECISION ONE, TWO
67: PARAMETER ( ONE = 1.0D+0, TWO = 2.0D+0 )
68: * ..
69: * .. Local Scalars ..
70: INTEGER I, J
71: DOUBLE PRECISION BIGNUM, EPS, SMLNUM, TEMP
72: * ..
73: * .. External Subroutines ..
74: EXTERNAL DLASWP, DSCAL
75: * ..
76: * .. External Functions ..
77: INTEGER IDAMAX
78: DOUBLE PRECISION DLAMCH
79: EXTERNAL IDAMAX, DLAMCH
80: * ..
81: * .. Intrinsic Functions ..
82: INTRINSIC ABS
83: * ..
84: * .. Executable Statements ..
85: *
86: * Set constant to control owerflow
87: *
88: EPS = DLAMCH( 'P' )
89: SMLNUM = DLAMCH( 'S' ) / EPS
90: BIGNUM = ONE / SMLNUM
91: CALL DLABAD( SMLNUM, BIGNUM )
92: *
93: * Apply permutations IPIV to RHS
94: *
95: CALL DLASWP( 1, RHS, LDA, 1, N-1, IPIV, 1 )
96: *
97: * Solve for L part
98: *
99: DO 20 I = 1, N - 1
100: DO 10 J = I + 1, N
101: RHS( J ) = RHS( J ) - A( J, I )*RHS( I )
102: 10 CONTINUE
103: 20 CONTINUE
104: *
105: * Solve for U part
106: *
107: SCALE = ONE
108: *
109: * Check for scaling
110: *
111: I = IDAMAX( N, RHS, 1 )
112: IF( TWO*SMLNUM*ABS( RHS( I ) ).GT.ABS( A( N, N ) ) ) THEN
113: TEMP = ( ONE / TWO ) / ABS( RHS( I ) )
114: CALL DSCAL( N, TEMP, RHS( 1 ), 1 )
115: SCALE = SCALE*TEMP
116: END IF
117: *
118: DO 40 I = N, 1, -1
119: TEMP = ONE / A( I, I )
120: RHS( I ) = RHS( I )*TEMP
121: DO 30 J = I + 1, N
122: RHS( I ) = RHS( I ) - RHS( J )*( A( I, J )*TEMP )
123: 30 CONTINUE
124: 40 CONTINUE
125: *
126: * Apply permutations JPIV to the solution (RHS)
127: *
128: CALL DLASWP( 1, RHS, LDA, 1, N-1, JPIV, -1 )
129: RETURN
130: *
131: * End of DGESC2
132: *
133: END
CVSweb interface <joel.bertrand@systella.fr>