1 SUBROUTINE DGESC2( N, A, LDA, RHS, IPIV, JPIV, SCALE )
2 *
3 * -- LAPACK auxiliary routine (version 3.2.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 * June 2010
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
2 *
3 * -- LAPACK auxiliary routine (version 3.2.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 * June 2010
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