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
|
// Jacobian.hh -- generic data structures for the Jacobian matrix
// $Header$
//
// Jacobian -- ABC to describe Jacobian matrix
// Jacobian_sparsity -- ABC to accumulate Jacobian sparsity info
// Jacobian_matrix -- ABC to store/use Jacobian matrix
//
// decode_Jacobian_type - decode string into internal enum
// new_Jacobian - factory method
//
//
// prerequisites:
// "../patch/patch_system.hh"
//
//******************************************************************************
//******************************************************************************
//******************************************************************************
//
// These classes are used to store and manipulate Jacobian matrices
// for a patch system.
//
// The inheritance diagram for the Jacobian classes looks like this:
// Jacobian
// Jacobian_sparsity
// dense_Jacobian_sparsity #ifdef HAVE_DENSE_JACOBIAN
// row_sparse_Jacobian_sparsity #ifdef HAVE_ROW_SPARSE_JACOBIAN
// Jacobian_matrix
// dense_Jacobian_matrix #ifdef HAVE_DENSE_JACOBIAN
// row_sparse_Jacobian_matrix #ifdef HAVE_ROW_SPARSE_JACOBIAN
// where the most-derived classes are all conditional on the noted #ifdef
// symbols.
//
// The Jacobian_sparsity objects are used to accumulate information
// on the Jacobian's sparsity pattern; this is then used in constructing
// the corresponding Jacobian_matrix object (which actually stores
// the Jacobian, and contains member functions to solve elliptic systems
// using the Jacobian as the LHS matrix).
//
// We assume a "traversal" interface for computing Jacobians, defined
// by the Jacobian API. Because this API is shared by (all) the
// Jacobian_sparsity and Jacobian_matrix classes, the same traversal
// routines (prototyped to update Jacobian objects) can be used both
// to record sparsity patterns (ignoring the actual floating-point values)
// and to store the Jacobian matrix.
//
// Thus the typical code to construct and use a Jacobian matrix looks
// like this:
// // prototype
// void trverse_Jacobian(const patch_system& ps, Jacobian& Jac);
//
// Jacobian_sparsity& JS = new_Jacobian_sparsity(Jac_type, ps);
// traverse_Jacobian(ps, JS);
// Jacobian_matrix& Jac = new_Jacobian_matrix(Jac_type, ps, JS);
// traverse_Jacobian(ps, Jac);
// Jac.solve_linear_system(...);
//
//
// A row/column index of the Jacobian (denoted II/JJ) is a 0-origin grid
// point number within the patch system.
//
// Note that the APIs here implicitly assume there is only a *single* gridfn
// in the Jacobian computation. (If we had N gridfns for this, then the
// Jacobian would really be a block matrix with N*N blocks.) This isn't
// very general, but matches our present use in this thorn.
//
//******************************************************************************
// ABC to define Jacobian-traversal interface
class Jacobian
{
public:
// basic meta-info
patch_system& my_patch_system() const { return ps_; }
int NN() const { return NN_; }
// convert (patch,irho,isigma) <--> row/column index
int II_of_patch_irho_isigma(const patch& p, int irho, int isigma)
const
{ return ps_.gpn_of_patch_irho_isigma(p, irho,isigma); }
const patch& patch_irho_isigma_of_II(int II, int& irho, int& isigma)
const
{ return ps_.patch_irho_isigma_of_gpn(II, irho,isigma); }
// zero the entire matrix
virtual void zero_matrix() = 0;
// set a matrix element to a specified value
virtual void set_element(int II, int JJ, fp value) = 0;
// sum a value into a matrix element
virtual void sum_into_element(int II, int JJ, fp value) = 0;
protected:
// constructor, destructor
Jacobian(patch_system& ps)
: ps_(ps),
NN_(ps.N_grid_points())
{ }
public:
virtual ~Jacobian() { }
private:
// we forbid copying and passing by value
// by declaring the copy constructor and assignment operator
// private, but never defining them
Jacobian(const Jacobian &rhs);
Jacobian& operator=(const Jacobian &rhs);
protected:
patch_system& ps_;
int NN_; // ps_.N_grid_points()
};
//******************************************************************************
// ABC to accumulate Jacobian sparsity info
class Jacobian_sparsity
: public Jacobian
{
public:
Jacobian_sparsity(patch_system& ps,
bool print_msg_flag = false)
: Jacobian(ps)
{ }
virtual ~Jacobian_sparsity() { }
};
//******************************************************************************
// Jacobian_matrix -- ABC to store/use Jacobian matrix
class Jacobian_matrix
: public Jacobian
{
public:
// access (get) a matrix element
virtual fp element(int II, int JJ) const = 0;
// is a given element explicitly stored, or implicitly 0 via sparsity
virtual bool is_explicitly_stored(int II, int JJ) const = 0;
// solve linear system J.x = rhs
// ... rhs and x are nominal-grid gridfns
// ... may modify Jacobian matrix (eg for LU decomposition)
// ... returns 0.0 if matrix is numerically singular
// condition number (> 0.0) if known
// -1.0 if condition number is unknown
virtual fp solve_linear_system(int rhs_gfn, int x_gfn) = 0;
// constructor, destructor
// ... note Jacobian_sparsity argument of constructor is non-const
// to allow us to take over ownership of data structures
Jacobian_matrix(patch_system& ps,
Jacobian_sparsity& JS,
bool print_msg_flag = false)
: Jacobian(ps)
{ }
virtual ~Jacobian_matrix() { }
};
//******************************************************************************
//******************************************************************************
//******************************************************************************
enum Jacobian_type
{
#ifdef HAVE_DENSE_JACOBIAN
Jacobian_type__dense_matrix,
#endif
#ifdef HAVE_ROW_SPARSE_JACOBIAN
Jacobian_type__row_sparse_matrix // no comma on last entry in enum
#endif
};
//
// prototypes of various Jacobian-related functions
//
// decode string into internal enum
enum Jacobian_type
decode_Jacobian_type(const char Jacobian_type_string[]);
// "object factory" routines to construct and return
// pointers to new-allocated objects of specified derived types
class Jacobian_sparsity;
class Jacobian_matrix;
Jacobian_sparsity* new_Jacobian_sparsity(enum Jacobian_type Jac_type,
patch_system& ps,
bool print_msg_flag = false);
Jacobian_matrix* new_Jacobian_matrix(enum Jacobian_type Jac_type,
patch_system& ps,
Jacobian_sparsity& JS,
bool print_msg_flag = false);
|
