aboutsummaryrefslogtreecommitdiff
path: root/src/elliptic/dense_Jacobian.hh
blob: 1e269d500a73200b8925eafcc8741b3efcad33d8 (plain)
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