1: *> \brief \b DGESC2
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 owerflow 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: *> \date November 2011
105: *
106: *> \ingroup doubleGEauxiliary
107: *
108: *> \par Contributors:
109: * ==================
110: *>
111: *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
112: *> Umea University, S-901 87 Umea, Sweden.
113: *
114: * =====================================================================
115: SUBROUTINE DGESC2( N, A, LDA, RHS, IPIV, JPIV, SCALE )
116: *
117: * -- LAPACK auxiliary routine (version 3.4.0) --
118: * -- LAPACK is a software package provided by Univ. of Tennessee, --
119: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
120: * November 2011
121: *
122: * .. Scalar Arguments ..
123: INTEGER LDA, N
124: DOUBLE PRECISION SCALE
125: * ..
126: * .. Array Arguments ..
127: INTEGER IPIV( * ), JPIV( * )
128: DOUBLE PRECISION A( LDA, * ), RHS( * )
129: * ..
130: *
131: * =====================================================================
132: *
133: * .. Parameters ..
134: DOUBLE PRECISION ONE, TWO
135: PARAMETER ( ONE = 1.0D+0, TWO = 2.0D+0 )
136: * ..
137: * .. Local Scalars ..
138: INTEGER I, J
139: DOUBLE PRECISION BIGNUM, EPS, SMLNUM, TEMP
140: * ..
141: * .. External Subroutines ..
142: EXTERNAL DLASWP, DSCAL
143: * ..
144: * .. External Functions ..
145: INTEGER IDAMAX
146: DOUBLE PRECISION DLAMCH
147: EXTERNAL IDAMAX, DLAMCH
148: * ..
149: * .. Intrinsic Functions ..
150: INTRINSIC ABS
151: * ..
152: * .. Executable Statements ..
153: *
154: * Set constant to control owerflow
155: *
156: EPS = DLAMCH( 'P' )
157: SMLNUM = DLAMCH( 'S' ) / EPS
158: BIGNUM = ONE / SMLNUM
159: CALL DLABAD( SMLNUM, BIGNUM )
160: *
161: * Apply permutations IPIV to RHS
162: *
163: CALL DLASWP( 1, RHS, LDA, 1, N-1, IPIV, 1 )
164: *
165: * Solve for L part
166: *
167: DO 20 I = 1, N - 1
168: DO 10 J = I + 1, N
169: RHS( J ) = RHS( J ) - A( J, I )*RHS( I )
170: 10 CONTINUE
171: 20 CONTINUE
172: *
173: * Solve for U part
174: *
175: SCALE = ONE
176: *
177: * Check for scaling
178: *
179: I = IDAMAX( N, RHS, 1 )
180: IF( TWO*SMLNUM*ABS( RHS( I ) ).GT.ABS( A( N, N ) ) ) THEN
181: TEMP = ( ONE / TWO ) / ABS( RHS( I ) )
182: CALL DSCAL( N, TEMP, RHS( 1 ), 1 )
183: SCALE = SCALE*TEMP
184: END IF
185: *
186: DO 40 I = N, 1, -1
187: TEMP = ONE / A( I, I )
188: RHS( I ) = RHS( I )*TEMP
189: DO 30 J = I + 1, N
190: RHS( I ) = RHS( I ) - RHS( J )*( A( I, J )*TEMP )
191: 30 CONTINUE
192: 40 CONTINUE
193: *
194: * Apply permutations JPIV to the solution (RHS)
195: *
196: CALL DLASWP( 1, RHS, LDA, 1, N-1, JPIV, -1 )
197: RETURN
198: *
199: * End of DGESC2
200: *
201: END
CVSweb interface <joel.bertrand@systella.fr>