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
|
// Jacobian.hh -- data structures for the Jacobian matrix
// $Header$
//
#ifdef HAVE_DENSE_JACOBIAN
// dense_Jacobian_sparsity -- (no-op) accumulate dense Jacobian sparsity info
// dense_Jacobian_matrix -- Jacobian stored as a dense matrix
#endif
//
//
// prerequisites:
// "../patch/patch_system.hh"
// "Jacobian.hh"
//
//******************************************************************************
#ifdef HAVE_DENSE_JACOBIAN
//
// This class accumulates the sparsity pattern of a dense-matrix Jacobian.
//
// Since a dense-matrix Jacobian doesn't care about the sparsity pattern,
// the member functions of this class are just no-ops.
//
class dense_Jacobian_sparsity
: public Jacobian_sparsity
{
public:
// record the zeroing of the entire matrix
void zero_matrix() { }
// record the setting of a matrix element to a specified value
void set_element(int II, int JJ, fp value) { }
// record the summing of a value into a matrix element
void sum_into_element(int II, int JJ, fp value) { }
// constructor, destructor
dense_Jacobian_sparsity(patch_system& ps,
bool print_msg_flag = false)
: Jacobian_sparsity(ps, print_msg_flag)
{ }
~dense_Jacobian_sparsity() { }
};
#endif // HAVE_DENSE_JACOBIAN
//******************************************************************************
#ifdef HAVE_DENSE_JACOBIAN
//
// This class stores the Jacobian as a dense matrix in Fortran (column)
// order.
//
class dense_Jacobian_matrix
: public Jacobian_matrix
{
public:
// zero the entire matrix
void zero_matrix();
// set a matrix element to a specified value
void set_element(int II, int JJ, fp value)
{ matrix_(JJ,II) = value; }
// sum a value into a matrix element
void sum_into_element(int II, int JJ, fp value)
{ matrix_(JJ,II) += value; }
// access (get) a matrix element
fp element(int II, int JJ)
const
{ return matrix_(JJ,II); }
// this is a dense matrix ==> all values are stored explicitly
bool is_explicitly_stored(int II, int JJ)
const
{ return true; }
// solve linear system J.x = rhs via LAPACK
// ... 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);
// constructor, destructor
public:
dense_Jacobian_matrix(patch_system& ps,
dense_Jacobian_sparsity& JS,
bool print_msg_flag = false);
~dense_Jacobian_matrix();
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 // HAVE_DENSE_JACOBIAN
|