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
|
// 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"
//******************************************************************************
//******************************************************************************
//******************************************************************************
//
// 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 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
default:
error_exit(ERROR_EXIT,
"new_Jacobian(): unknown method=(int)%d!\n",
int(Jac_method)); /*NOTREACHED*/
}
}
|
