1: *> \brief \b DGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed by sgetc2.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DGESC2 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesc2.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesc2.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesc2.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DGESC2( N, A, LDA, RHS, IPIV, JPIV, SCALE )
22: *
23: * .. Scalar Arguments ..
24: * INTEGER LDA, N
25: * DOUBLE PRECISION SCALE
26: * ..
27: * .. Array Arguments ..
28: * INTEGER IPIV( * ), JPIV( * )
29: * DOUBLE PRECISION A( LDA, * ), RHS( * )
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> DGESC2 solves a system of linear equations
39: *>
40: *> A * X = scale* RHS
41: *>
42: *> with a general N-by-N matrix A using the LU factorization with
43: *> complete pivoting computed by DGETC2.
44: *> \endverbatim
45: *
46: * Arguments:
47: * ==========
48: *
49: *> \param[in] N
50: *> \verbatim
51: *> N is INTEGER
52: *> The order of the matrix A.
53: *> \endverbatim
54: *>
55: *> \param[in] A
56: *> \verbatim
57: *> A is DOUBLE PRECISION array, dimension (LDA,N)
58: *> On entry, the LU part of the factorization of the n-by-n
59: *> matrix A computed by DGETC2: A = P * L * U * Q
60: *> \endverbatim
61: *>
62: *> \param[in] LDA
63: *> \verbatim
64: *> LDA is INTEGER
65: *> The leading dimension of the array A. LDA >= max(1, N).
66: *> \endverbatim
67: *>
68: *> \param[in,out] RHS
69: *> \verbatim
70: *> RHS is DOUBLE PRECISION array, dimension (N).
71: *> On entry, the right hand side vector b.
72: *> On exit, the solution vector X.
73: *> \endverbatim
74: *>
75: *> \param[in] IPIV
76: *> \verbatim
77: *> IPIV is INTEGER array, dimension (N).
78: *> The pivot indices; for 1 <= i <= N, row i of the
79: *> matrix has been interchanged with row IPIV(i).
80: *> \endverbatim
81: *>
82: *> \param[in] JPIV
83: *> \verbatim
84: *> JPIV is INTEGER array, dimension (N).
85: *> The pivot indices; for 1 <= j <= N, column j of the
86: *> matrix has been interchanged with column JPIV(j).
87: *> \endverbatim
88: *>
89: *> \param[out] SCALE
90: *> \verbatim
91: *> SCALE is DOUBLE PRECISION
92: *> On exit, SCALE contains the scale factor. SCALE is chosen
93: *> 0 <= SCALE <= 1 to prevent overflow in the solution.
94: *> \endverbatim
95: *
96: * Authors:
97: * ========
98: *
99: *> \author Univ. of Tennessee
100: *> \author Univ. of California Berkeley
101: *> \author Univ. of Colorado Denver
102: *> \author NAG Ltd.
103: *
104: *> \ingroup doubleGEauxiliary
105: *
106: *> \par Contributors:
107: * ==================
108: *>
109: *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
110: *> Umea University, S-901 87 Umea, Sweden.
111: *
112: * =====================================================================
113: SUBROUTINE DGESC2( N, A, LDA, RHS, IPIV, JPIV, SCALE )
114: *
115: * -- LAPACK auxiliary routine --
116: * -- LAPACK is a software package provided by Univ. of Tennessee, --
117: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
118: *
119: * .. Scalar Arguments ..
120: INTEGER LDA, N
121: DOUBLE PRECISION SCALE
122: * ..
123: * .. Array Arguments ..
124: INTEGER IPIV( * ), JPIV( * )
125: DOUBLE PRECISION A( LDA, * ), RHS( * )
126: * ..
127: *
128: * =====================================================================
129: *
130: * .. Parameters ..
131: DOUBLE PRECISION ONE, TWO
132: PARAMETER ( ONE = 1.0D+0, TWO = 2.0D+0 )
133: * ..
134: * .. Local Scalars ..
135: INTEGER I, J
136: DOUBLE PRECISION BIGNUM, EPS, SMLNUM, TEMP
137: * ..
138: * .. External Subroutines ..
139: EXTERNAL DLASWP, DSCAL, DLABAD
140: * ..
141: * .. External Functions ..
142: INTEGER IDAMAX
143: DOUBLE PRECISION DLAMCH
144: EXTERNAL IDAMAX, DLAMCH
145: * ..
146: * .. Intrinsic Functions ..
147: INTRINSIC ABS
148: * ..
149: * .. Executable Statements ..
150: *
151: * Set constant to control overflow
152: *
153: EPS = DLAMCH( 'P' )
154: SMLNUM = DLAMCH( 'S' ) / EPS
155: BIGNUM = ONE / SMLNUM
156: CALL DLABAD( SMLNUM, BIGNUM )
157: *
158: * Apply permutations IPIV to RHS
159: *
160: CALL DLASWP( 1, RHS, LDA, 1, N-1, IPIV, 1 )
161: *
162: * Solve for L part
163: *
164: DO 20 I = 1, N - 1
165: DO 10 J = I + 1, N
166: RHS( J ) = RHS( J ) - A( J, I )*RHS( I )
167: 10 CONTINUE
168: 20 CONTINUE
169: *
170: * Solve for U part
171: *
172: SCALE = ONE
173: *
174: * Check for scaling
175: *
176: I = IDAMAX( N, RHS, 1 )
177: IF( TWO*SMLNUM*ABS( RHS( I ) ).GT.ABS( A( N, N ) ) ) THEN
178: TEMP = ( ONE / TWO ) / ABS( RHS( I ) )
179: CALL DSCAL( N, TEMP, RHS( 1 ), 1 )
180: SCALE = SCALE*TEMP
181: END IF
182: *
183: DO 40 I = N, 1, -1
184: TEMP = ONE / A( I, I )
185: RHS( I ) = RHS( I )*TEMP
186: DO 30 J = I + 1, N
187: RHS( I ) = RHS( I ) - RHS( J )*( A( I, J )*TEMP )
188: 30 CONTINUE
189: 40 CONTINUE
190: *
191: * Apply permutations JPIV to the solution (RHS)
192: *
193: CALL DLASWP( 1, RHS, LDA, 1, N-1, JPIV, -1 )
194: RETURN
195: *
196: * End of DGESC2
197: *
198: END
CVSweb interface <joel.bertrand@systella.fr>