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
|
// Jacobian.cc -- generic routines for the Jacobian matrix
// $Header$
//
// <<<literal contents of "lapack.hh">>>
//
// decode_Jacobian_store_solve_method -- decode string into internal enum
// new_Jacobian -- object factory for Jacobian objects
//
#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"
// all the code in this file is inside this namespace
namespace AHFinderDirect
{
//******************************************************************************
//******************************************************************************
//******************************************************************************
//
// This function decodes a character string specifying a specific type
// (derived class) of Jacobian, into an internal enum.
//
enum Jacobian_store_solve_method
decode_Jacobian_store_solve_method
(const char Jacobian_store_solve_method_string[])
{
if (STRING_EQUAL(Jacobian_store_solve_method_string,
"dense matrix/LAPACK"))
then {
#ifdef HAVE_DENSE_JACOBIAN__LAPACK
return Jacobian__dense_matrix__LAPACK;
#endif
}
else if (STRING_EQUAL(Jacobian_store_solve_method_string,
"row-oriented sparse matrix/ILUCG"))
then {
#ifdef HAVE_ROW_SPARSE_JACOBIAN__ILUCG
return Jacobian__row_sparse_matrix__ILUCG;
#endif
}
else if (STRING_EQUAL(Jacobian_store_solve_method_string,
"row-oriented sparse matrix/UMFPACK"))
then {
#ifdef HAVE_ROW_SPARSE_JACOBIAN__UMFPACK
return Jacobian__row_sparse_matrix__UMFPACK;
#endif
}
else error_exit(ERROR_EXIT,
"decode_Jacobian_store_solve_method():\n"
" unknown Jacobian_store_solve_method_string=\"%s\"!\n",
Jacobian_store_solve_method_string); /*NOTREACHED*/
// fall through to here ==> we recognize the matrix store_solve_method,
// but it's not configured
error_exit(ERROR_EXIT,
"\n"
" decode_Jacobian_store_solve_method():\n"
" Jacobian store_solve_method=\"%s\"\n"
" is not configured in this binary; see \"src/include/config.hh\"\n"
" for details on what methods are configured and how to change this\n"
,
Jacobian_store_solve_method_string); /*NOTREACHED*/
}
//******************************************************************************
//
// This function is an "object factory" for Jacobian objects: it constructs
// and returns a pointer to a new-allocated Jacobian 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* new_Jacobian(enum Jacobian_store_solve_method Jac_method,
patch_system& ps,
bool print_msg_flag /* = false */)
{
switch (Jac_method)
{
#ifdef HAVE_DENSE_JACOBIAN__LAPACK
case Jacobian__dense_matrix__LAPACK:
return new dense_Jacobian__LAPACK(ps, print_msg_flag);
#endif
#ifdef HAVE_ROW_SPARSE_JACOBIAN__ILUCG
case Jacobian__row_sparse_matrix__ILUCG:
return new row_sparse_Jacobian__ILUCG(ps, print_msg_flag);
#endif
#ifdef HAVE_ROW_SPARSE_JACOBIAN__UMFPACK
case Jacobian__row_sparse_matrix__UMFPACK:
return new row_sparse_Jacobian__UMFPACK(ps, print_msg_flag);
#endif
default:
error_exit(ERROR_EXIT,
"new_Jacobian(): unknown method=(int)%d!\n",
int(Jac_method)); /*NOTREACHED*/
}
}
//******************************************************************************
//******************************************************************************
//******************************************************************************
} // namespace AHFinderDirect
|
