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
|
// Jacobian.hh -- data structures for the Jacobian matrix
// $Header$
//
// Jacobian -- abstract base class to describe a Jacobian matrix
#ifdef HAVE_DENSE_JACOBIAN
// dense_Jacobian -- Jacobian stored as a dense matrix
#endif
//
// decode_Jacobian_type - decode string into internal enum
// new_Jacobian - factory method
//
//
// prerequisites:
// "stdc.h"
// "config.hh"
// "../jtutil/util.hh"
// "../jtutil/array.hh"
// "../jtutil/cpm_map.hh"
// "../jtutil/linear_map.hh"
// "coords.hh"
// "grid.hh"
// "fd_grid.hh"
// "patch.hh"
// "patch_edge.hh"
// "patch_interp.hh"
// "ghost_zone.hh"
// "patch_system.hh"
//
//******************************************************************************
//
// A Jacobian object stores the (a) Jacobian matrix for a patch system.
// Jacobian is an abstract base class, from which we derive a separate
// concrete class for each type of sparsity pattern (more precisely for
// each type of storage scheme) of the Jacobian matrix.
//
// 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.
//
class Jacobian
{
public:
// access a matrix element via row/column indices
virtual const fp& operator()(int II, int JJ) const = 0; // rvalue
virtual fp& operator()(int II, int JJ) = 0; // lvalue
// access a matrix element via row index and column (patch,irho,isigma)
const fp& operator() // rvalue
(int II, const patch& JJ_patch, int JJ_irho, int JJ_isigma)
const
{
const int JJ
= ps_.gpn_of_patch_irho_isigma(JJ_patch, JJ_irho,JJ_isigma);
return operator()(II,JJ);
}
fp& operator() // lvalue
(int II, const patch& JJ_patch, int JJ_irho, int JJ_isigma)
{
const int JJ
= ps_.gpn_of_patch_irho_isigma(JJ_patch, JJ_irho,JJ_isigma);
return operator()(II,JJ);
}
// is a given element explicitly stored, or implicitly 0 via sparsity
virtual bool is_explicitly_stored(int II, int JJ) const = 0;
// zero the whole matrix or a single row
virtual void zero_matrix() = 0;
virtual void zero_row(int II) = 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;
// basic meta-info
patch_system& my_patch_system() const { return ps_; }
int NN() const { return NN_; }
// constructor, destructor
// ... for analyze-factor-solve sparse-matrix derived classes,
// construction may include the "analyze" operation
Jacobian(patch_system& ps);
virtual ~Jacobian() { }
protected:
patch_system& ps_;
int NN_; // ps_.N_grid_points()
};
//******************************************************************************
#ifdef HAVE_DENSE_JACOBIAN
//
// This class stores the Jacobian as a dense matrix in Fortran (column)
// order.
//
class dense_Jacobian
: public Jacobian
{
public:
// access a matrix element via row/column indices
const fp& operator()(int II, int JJ) const // rvalue
{ return matrix_(JJ,II); }
fp& operator()(int II, int JJ) // lvalue
{ return matrix_(JJ,II); }
bool is_explicitly_stored(int II, int JJ) const
{ return true; }
// zero the whole matrix or a single row
void zero_matrix();
void zero_row(int II);
// solve linear system J.x = rhs via LAPACK
// ... rhs and x are nominal-grid gridfns
// ... overwrites Jacobian matrix with its LU decomposition
// ... returns 0.0 if matrix is numerically singular, or
// condition number (> 0.0)
fp solve_linear_system(int rhs_gfn, int x_gfn);
// constructor, destructor
public:
dense_Jacobian(patch_system& ps);
~dense_Jacobian();
private:
// subscripts are (JJ,II) (both 0-origin)
jtutil::array2d<fp> matrix_;
// pivot vector for LAPACK routines
integer *pivot_;
// work vectors for LAPACK condition number computation
integer *iwork_;
fp *rwork_;
};
#endif
//******************************************************************************
enum Jacobian_type
{
#ifdef HAVE_DENSE_JACOBIAN
Jacobian_type__dense_matrix // no comma on last entry in enum
#endif
};
// decode string into internal enum
enum Jacobian_type
decode_Jacobian_type(const char Jacobian_type_string[]);
// construct and return new-allocated Jacobian object of specified type
Jacobian& new_Jacobian(patch_system& ps, enum Jacobian_type Jac_type);
|
