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
|
// Jacobian.cc -- generic routines for the Jacobian matrix
// $Header$
//
// <<<literal contents of "lapack.hh">>>
//
// decode_Jacobian_type
// new_Jacobian_sparsity
// new_Jacobian_matrix
//
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include <string.h>
#include "util_Table.h"
#include "cctk.h"
#include "cctk_Arguments.h"
#include "config.h"
#include "stdc.h"
#include "../jtutil/util.hh"
#include "../jtutil/array.hh"
#include "../jtutil/cpm_map.hh"
#include "../jtutil/linear_map.hh"
using jtutil::error_exit;
#include "../patch/coords.hh"
#include "../patch/grid.hh"
#include "../patch/fd_grid.hh"
#include "../patch/patch.hh"
#include "../patch/patch_edge.hh"
#include "../patch/patch_interp.hh"
#include "../patch/ghost_zone.hh"
#include "../patch/patch_system.hh"
#include "Jacobian.hh"
#include "dense_Jacobian.hh"
#include "row_sparse_Jacobian.hh"
//******************************************************************************
//******************************************************************************
//******************************************************************************
//
// This function decodes a character string specifying a type (derived class)
// of Jacobian, into an internal enum.
//
enum Jacobian_type
decode_Jacobian_type(const char Jacobian_type_string[])
{
if (STRING_EQUAL(Jacobian_type_string, "dense matrix"))
then {
#ifdef HAVE_DENSE_JACOBIAN
return Jacobian_type__dense_matrix;
#endif
}
else if (STRING_EQUAL(Jacobian_type_string, "row-oriented sparse matrix"))
then {
#ifdef HAVE_ROW_SPARSE_JACOBIAN
return Jacobian_type__row_sparse_matrix;
#endif
}
else error_exit(ERROR_EXIT,
"decode_Jacobian_type(): unknown Jacobian_type_string=\"%s\"!\n",
Jacobian_type_string); /*NOTREACHED*/
// fall through to here ==> we recognize the matrix type,
// but it's not configured
error_exit(ERROR_EXIT,
"\n"
" decode_Jacobian_type():\n"
" Jacobian type \"%s\" is not configured in this binary;\n"
" see \"src/include/config.hh\" for details on what Jacobian types\n"
" are configured and how to change this\n"
,
Jacobian_type_string); /*NOTREACHED*/
}
//******************************************************************************
//
// This function is an "object factory" for Jacobian_sparsity objects:
// it constructs and returns a pointer to a new-allocated Jacobian_sparsity
// object of the specified derived type.
//
// FIXME: the patch system shouldn't really have to be non-const, but
// the Jacobian constructors all require this to allow the
// linear solvers to directly update gridfns
//
Jacobian_sparsity* new_Jacobian_sparsity(enum Jacobian_type Jac_type,
patch_system& ps,
bool print_msg_flag /* = false */)
{
switch (Jac_type)
{
#ifdef HAVE_DENSE_JACOBIAN
case Jacobian_type__dense_matrix:
return new dense_Jacobian_sparsity(ps, print_msg_flag);
#endif
#ifdef HAVE_ROW_SPARSE_JACOBIAN
case Jacobian_type__row_sparse_matrix:
return new row_sparse_Jacobian_sparsity(ps, print_msg_flag);
#endif
default:
error_exit(ERROR_EXIT,
"new_Jacobian_sparsity(): unknown Jacobian_type=(int)%d!\n",
Jac_type); /*NOTREACHED*/
}
}
//******************************************************************************
//
// This function is an "object factory" for Jacobian_matrix objects:
// it constructs and returns a pointer to a new-allocated Jacobian_matrix
// object of the specified derived type.
//
// FIXME: the patch system shouldn't really have to be non-const, but
// the Jacobian constructors all require this to allow the
// linear solvers to directly update gridfns
//
Jacobian_matrix* new_Jacobian_matrix(enum Jacobian_type Jac_type,
patch_system& ps,
Jacobian_sparsity& JS,
bool print_msg_flag /* = false */)
{
switch (Jac_type)
{
#ifdef HAVE_DENSE_JACOBIAN
case Jacobian_type__dense_matrix:
return new dense_Jacobian_matrix
(ps, static_cast<dense_Jacobian_sparsity&>(JS),
print_msg_flag);
#endif
#ifdef HAVE_ROW_SPARSE_JACOBIAN
case Jacobian_type__row_sparse_matrix:
return new row_sparse_Jacobian_matrix
(ps, static_cast<row_sparse_Jacobian_sparsity&>(JS),
print_msg_flag);
#endif
default:
error_exit(ERROR_EXIT,
"new_Jacobian_matrix(): unknown Jacobian_type=(int)%d!\n",
Jac_type); /*NOTREACHED*/
}
}
|
