Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
ElmerCSC
GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/mathlibs/src/lapack/cgesc2.f
5194 views
1
SUBROUTINE CGESC2( N, A, LDA, RHS, IPIV, JPIV, SCALE )
2
*
3
* -- LAPACK auxiliary routine (version 3.0) --
4
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5
* Courant Institute, Argonne National Lab, and Rice University
6
* June 30, 1999
7
*
8
* .. Scalar Arguments ..
9
INTEGER LDA, N
10
REAL SCALE
11
* ..
12
* .. Array Arguments ..
13
INTEGER IPIV( * ), JPIV( * )
14
COMPLEX A( LDA, * ), RHS( * )
15
* ..
16
*
17
* Purpose
18
* =======
19
*
20
* CGESC2 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 CGETC2.
26
*
27
*
28
* Arguments
29
* =========
30
*
31
* N (input) INTEGER
32
* The number of columns of the matrix A.
33
*
34
* A (input) COMPLEX array, dimension (LDA, N)
35
* On entry, the LU part of the factorization of the n-by-n
36
* matrix A computed by CGETC2: A = P * L * U * Q
37
*
38
* LDA (input) INTEGER
39
* The leading dimension of the array A. LDA >= max(1, N).
40
*
41
* RHS (input/output) COMPLEX array, dimension N.
42
* On entry, the right hand side vector b.
43
* On exit, the solution vector X.
44
*
45
* IPIV (iput) INTEGER array, dimension (N).
46
* The pivot indices; for 1 <= i <= N, row i of the
47
* matrix has been interchanged with row IPIV(i).
48
*
49
* JPIV (iput) INTEGER array, dimension (N).
50
* The pivot indices; for 1 <= j <= N, column j of the
51
* matrix has been interchanged with column JPIV(j).
52
*
53
* SCALE (output) REAL
54
* On exit, SCALE contains the scale factor. SCALE is chosen
55
* 0 <= SCALE <= 1 to prevent owerflow in the solution.
56
*
57
* Further Details
58
* ===============
59
*
60
* Based on contributions by
61
* Bo Kagstrom and Peter Poromaa, Department of Computing Science,
62
* Umea University, S-901 87 Umea, Sweden.
63
*
64
* =====================================================================
65
*
66
* .. Parameters ..
67
REAL ZERO, ONE, TWO
68
PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 )
69
* ..
70
* .. Local Scalars ..
71
INTEGER I, J
72
REAL BIGNUM, EPS, SMLNUM
73
COMPLEX TEMP
74
* ..
75
* .. External Subroutines ..
76
EXTERNAL CLASWP, CSCAL, SLABAD
77
* ..
78
* .. External Functions ..
79
INTEGER ICAMAX
80
REAL SLAMCH
81
EXTERNAL ICAMAX, SLAMCH
82
* ..
83
* .. Intrinsic Functions ..
84
INTRINSIC ABS, CMPLX, REAL
85
* ..
86
* .. Executable Statements ..
87
*
88
* Set constant to control overflow
89
*
90
EPS = SLAMCH( 'P' )
91
SMLNUM = SLAMCH( 'S' ) / EPS
92
BIGNUM = ONE / SMLNUM
93
CALL SLABAD( SMLNUM, BIGNUM )
94
*
95
* Apply permutations IPIV to RHS
96
*
97
CALL CLASWP( 1, RHS, LDA, 1, N-1, IPIV, 1 )
98
*
99
* Solve for L part
100
*
101
DO 20 I = 1, N - 1
102
DO 10 J = I + 1, N
103
RHS( J ) = RHS( J ) - A( J, I )*RHS( I )
104
10 CONTINUE
105
20 CONTINUE
106
*
107
* Solve for U part
108
*
109
SCALE = ONE
110
*
111
* Check for scaling
112
*
113
I = ICAMAX( N, RHS, 1 )
114
IF( TWO*SMLNUM*ABS( RHS( I ) ).GT.ABS( A( N, N ) ) ) THEN
115
TEMP = CMPLX( ONE / TWO, ZERO ) / ABS( RHS( I ) )
116
CALL CSCAL( N, TEMP, RHS( 1 ), 1 )
117
SCALE = SCALE*REAL( TEMP )
118
END IF
119
DO 40 I = N, 1, -1
120
TEMP = CMPLX( ONE, ZERO ) / A( I, I )
121
RHS( I ) = RHS( I )*TEMP
122
DO 30 J = I + 1, N
123
RHS( I ) = RHS( I ) - RHS( J )*( A( I, J )*TEMP )
124
30 CONTINUE
125
40 CONTINUE
126
*
127
* Apply permutations JPIV to the solution (RHS)
128
*
129
CALL CLASWP( 1, RHS, LDA, 1, N-1, JPIV, -1 )
130
RETURN
131
*
132
* End of CGESC2
133
*
134
END
135
136