1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
|
// Jacobian.hh -- data structures for the Jacobian matrix
// $Header$
//
#ifdef HAVE_ROW_SPARSE_JACOBIAN
// row_sparse_Jacobian -- Jacobian stored as row-oriented sparse matrix
#endif
#ifdef HAVE_ROW_SPARSE_JACOBIAN__ILUCG
// row_sparse_Jacobian__ILUCG -- row_sparse_Jacobian with ILUCG linear solver
#endif
//
#ifndef AHFINDERDIRECT__ROW_SPARSE_JACOBIAN_HH
#define AHFINDERDIRECT__ROW_SPARSE_JACOBIAN_HH
//
// prerequisites:
// "../patch/patch_system.hh"
// "Jacobian.hh"
//
#undef DEBUG_ROW_SPARSE_JACOBIAN // defined this for *very* detailed
// printing of data structures
//******************************************************************************
#ifdef HAVE_ROW_SPARSE_JACOBIAN
//
// This class stores the Jacobian as a row-oriented sparse matrix.
//
// This sparse matrix format is widely used by sparse matrix libraries.
// The format is as follows: The matrix is represented by 3 arrays
// (where () denote Fortran 1-origin subscripting):
// IA(N_rows+1) = Fortran (1-origin) indices in A of the start
// of each row's nonzeros; IA(N_rows+1) = N_nonzeros + 1
// JA(N_nonzeros) = Fortran (1-origin) column numbers of the
// corresponding entries in A
// A(N_nonzeros) = nonzero matrix elements, ordered by rows;
// within a row the matrix elements may be in
// any order
// IA and JA are arrays of integer so as to be storage-compatible with
// Fortran.
//
// For example, the matrix
// [ 11.0 12.0 13.0 16.0 ]
// [ 21.0 22.0 24.0 25.0 ]
// [ 32.0 34.0 ]
// [ 42.0 44.0 46.0 ]
// [ 53.0 55.0 ]
// [ 61.0 64.0 66.0 ]
// could be represented by the arrays (Fortran indexing)
// IA( 1) = 1
// JA( 1) = 1 JA( 2) = 2 JA( 3) = 3 JA( 4) = 6
// A( 1) = 11.0 A( 2) = 12.0 A( 3) = 13.0 A( 4) = 16.0
// IA( 2) = 5
// JA( 5) = 1 JA( 6) = 2 JA( 7) = 4 JA( 8) = 5
// A( 5) = 21.0 A( 6) = 22.0 A( 7) = 24.0 A( 8) = 25.0
// IA( 3) = 9
// JA( 9) = 2 JA(10) = 4
// A( 9) = 32.0 A(10) = 34.0
// IA( 4) = 11
// JA(11) = 2 JA(12) = 4 JA(13) = 6
// A(11) = 42.0 A(12) = 44.0 A(13) = 46.0
// IA( 5) = 14
// JA(14) = 3 JA(15) = 5
// A(14) = 53.0 A(15) = 55.0
// IA( 6) = 16
// JA(16) = 1 JA(17) = 4 JA(18) = 6
// A(16) = 61.0 A(17) = 64.0 A(18) = 66.0
// IA( 7) = 19
//
class row_sparse_Jacobian
: public Jacobian
{
public:
//
// routines to access the matrix
//
// get a matrix element
fp element(int II, int JJ) const;
// is the matrix element (II,JJ) stored explicitly?
bool is_explicitly_stored(int II, int JJ) const
{ return find_element(II,JJ) > 0; }
//
// routines for setting values in the matrix
//
// zero the entire matrix
void zero_matrix();
// set a matrix element to a specified value
void set_element(int II, int JJ, fp value);
// sum a value into a matrix element
void sum_into_element(int II, int JJ, fp value);
//
// constructor, destructor
//
protected:
// the constructor only uses ps to get the size of the matrix
// ... print_msg_flag controls messages about data structure setup
// during construction and element-setting process
row_sparse_Jacobian(patch_system& ps,
bool print_msg_flag = false);
~row_sparse_Jacobian();
//
// internal routines
//
protected:
// return position (1-origin) of (II,JJ) in A_ and JA_
// return 0 if not found
int find_element(int II, int JJ) const;
// insert new array element (II,JJ) in matrix,
// starting new row and/or growing arrays if necessary
// return position (1-origin) in A_ and JA_
int insert_element(int II, int JJ, fp value);
#ifdef DEBUG_ROW_SPARSE_JACOBIAN
// print data structure (IA_ and JA_ arrays)
void print_data_structure() const;
#endif
// access IA[] and JA[] and A[] using 1-origin Fortran indices
int& fIA(int fII) const
{
assert(fII >= 1);
assert(fII <= current_N_rows_+1);
return IA_[csub(fII)];
}
int& fJA(int fposn) const
{
assert(fposn >= 1);
assert(fposn <= N_nonzeros_);
return JA_[csub(fposn)];
}
fp& fA(int fposn) const
{
assert(fposn >= 1);
assert(fposn <= N_nonzeros_);
return A_[csub(fposn)];
}
protected:
// actual size, size allocated for JA_ and A_
int N_nonzeros_, N_nonzeros_allocated_;
// at present rows 1 <= fii <= current_N_rows_ are valid
int current_N_rows_;
// --> new[]-allocated arrays
// ***** due to constructor initialization list ordering,
// ***** these must be declared *AFTER* N_nonzeros_allocated_
integer* IA_; // size N_rows_+1
// valid for Fortran indices
// 1 <= fII <= current_N_rows_+1
// ... these are grown (reallocated) as necessary by insert_element()
// ... using STL vector would be much cleaner here, but alas
// there are too many broken compilers/linkers in the world
// world don't fully grok vector :( :(
integer* JA_; // size N_nonzeros_allocated_
// valid for Fortran indices
// 1 <= fposn < IA(current_N_rows_)
fp* A_; // ditto
};
#endif // HAVE_ROW_SPARSE_JACOBIAN
//******************************************************************************
#ifdef HAVE_ROW_SPARSE_JACOBIAN__ILUCG
//
// This class defines the linear solver routine using ILUCG.
//
class row_sparse_Jacobian__ILUCG
: public row_sparse_Jacobian
{
public:
// solve linear system J.x = rhs via ILUCG
// ... rhs and x are nominal-grid gridfns
// ... overwrites Jacobian matrix with its LU decomposition
// ... returns condition number if known,
// 0.0 if matrix is numerically singular,
// -1.0 if condition number is unknown
fp solve_linear_system(int rhs_gfn, int x_gfn,
const struct linear_solver_pars& pars,
bool print_msg_flag);
// constructor, destructor
public:
// the constructor only uses ps to get the size of the matrix
row_sparse_Jacobian__ILUCG(patch_system& ps,
bool print_msg_flag = false);
~row_sparse_Jacobian__ILUCG();
private:
bool print_msg_flag_;
// work vectors for ILUCG subroutine
// ... allocated by solve_linear_system() on first call
integer* itemp_;
fp* rtemp_;
};
#endif // HAVE_ROW_SPARSE_JACOBIAN__ILUCG
//******************************************************************************
#endif // AHFINDERDIRECT__ROW_SPARSE_JACOBIAN_HH
|
