-
-
Notifications
You must be signed in to change notification settings - Fork 0
/
octave-sundials6.patch
203 lines (185 loc) · 7.54 KB
/
octave-sundials6.patch
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
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
diff -up octave-6.4.0/libinterp/dldfcn/__ode15__.cc.1~ octave-6.4.0/libinterp/dldfcn/__ode15__.cc
--- octave-6.4.0/libinterp/dldfcn/__ode15__.cc.1~ 2021-10-30 16:20:24.000000000 +0200
+++ octave-6.4.0/libinterp/dldfcn/__ode15__.cc 2021-12-30 22:33:41.057334286 +0100
@@ -188,6 +188,9 @@ namespace octave
IDAFree (&m_mem);
SUNLinSolFree (m_sunLinearSolver);
SUNMatDestroy (m_sunJacMatrix);
+# if defined (HAVE_SUNDIALS_SUNCONTEXT)
+ SUNContext_Free (&m_sunContext);
+# endif
}
IDA&
@@ -247,7 +250,11 @@ namespace octave
static ColumnVector NVecToCol (N_Vector& v, long int n);
+# if defined (HAVE_SUNDIALS_SUNCONTEXT)
+ N_Vector ColToNVec (const ColumnVector& data, long int n);
+# else
static N_Vector ColToNVec (const ColumnVector& data, long int n);
+# endif
void
set_up (const ColumnVector& y);
@@ -356,6 +363,9 @@ namespace octave
DAEJacFuncSparse m_jacspfun;
DAEJacCellDense m_jacdcell;
DAEJacCellSparse m_jacspcell;
+# if defined (HAVE_SUNDIALS_SUNCONTEXT)
+ SUNContext m_sunContext;
+# endif
SUNMatrix m_sunJacMatrix;
SUNLinearSolver m_sunLinearSolver;
};
@@ -385,6 +395,12 @@ namespace octave
puntrr[i] = res(i);
}
+# if defined (HAVE_SUNDIALS_SUNCONTEXT)
+# define OCTAVE_SUNCONTEXT , m_sunContext
+# else
+# define OCTAVE_SUNCONTEXT
+# endif
+
void
IDA::set_up (const ColumnVector& y)
{
@@ -396,11 +412,11 @@ namespace octave
// FIXME : one should not allocate space for a full Jacobian
// when using a sparse format. Consider allocating less space
// then possibly using SUNSparseMatrixReallocate to increase it.
- m_sunJacMatrix = SUNSparseMatrix (m_num, m_num, m_num*m_num, CSC_MAT);
+ m_sunJacMatrix = SUNSparseMatrix (m_num, m_num, m_num*m_num, CSC_MAT OCTAVE_SUNCONTEXT);
if (! m_sunJacMatrix)
error ("Unable to create sparse Jacobian for Sundials");
- m_sunLinearSolver = SUNLinSol_KLU (yy, m_sunJacMatrix);
+ m_sunLinearSolver = SUNLinSol_KLU (yy, m_sunJacMatrix OCTAVE_SUNCONTEXT);
if (! m_sunLinearSolver)
error ("Unable to create KLU sparse solver");
@@ -417,11 +433,12 @@ namespace octave
else
{
- m_sunJacMatrix = SUNDenseMatrix (m_num, m_num);
+ m_sunJacMatrix = SUNDenseMatrix (m_num, m_num OCTAVE_SUNCONTEXT);
if (! m_sunJacMatrix)
error ("Unable to create dense Jacobian for Sundials");
- m_sunLinearSolver = SUNLinSol_Dense (yy, m_sunJacMatrix);
+ m_sunLinearSolver = SUNLinSol_Dense (yy, m_sunJacMatrix
+ OCTAVE_SUNCONTEXT);
if (! m_sunLinearSolver)
error ("Unable to create dense linear solver");
@@ -505,7 +522,7 @@ namespace octave
N_Vector
IDA::ColToNVec (const ColumnVector& data, long int n)
{
- N_Vector v = N_VNew_Serial (n);
+ N_Vector v = N_VNew_Serial (n OCTAVE_SUNCONTEXT);
realtype *punt = nv_data_s (v);
@@ -528,7 +545,13 @@ namespace octave
IDA::initialize (void)
{
m_num = m_y0.numel ();
+# if defined (HAVE_SUNDIALS_SUNCONTEXT)
+ if (SUNContext_Create (nullptr, &m_sunContext) < 0)
+ error ("__ode15__: unable to create context for SUNDIALS");
+ m_mem = IDACreate (m_sunContext);
+# else
m_mem = IDACreate ();
+# endif
N_Vector yy = ColToNVec (m_y0, m_num);
@@ -687,7 +710,7 @@ namespace octave
if (status == 0)
{
// Interpolate in tend
- N_Vector dky = N_VNew_Serial (m_num);
+ N_Vector dky = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT);
if (IDAGetDky (m_mem, tend, 0, dky) != 0)
error ("IDA failed to interpolate y");
@@ -837,9 +860,9 @@ namespace octave
realtype h = 0, tcur = 0;
bool status = 0;
- N_Vector dky = N_VNew_Serial (m_num);
+ N_Vector dky = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT);
- N_Vector dkyp = N_VNew_Serial (m_num);
+ N_Vector dkyp = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT);
ColumnVector yout (m_num);
ColumnVector ypout (m_num);
diff -up octave-6.4.0/m4/acinclude.m4.1~ octave-6.4.0/m4/acinclude.m4
--- octave-6.4.0/m4/acinclude.m4.1~ 2021-10-30 16:20:24.000000000 +0200
+++ octave-6.4.0/m4/acinclude.m4 2021-12-30 22:35:21.123832002 +0100
@@ -2696,6 +2696,7 @@ AC_DEFUN([OCTAVE_CHECK_SUNDIALS_COMPATIB
octave_have_sundials_compatible_api=no
fi
AC_MSG_RESULT([$octave_have_sundials_compatible_api])
+ AC_DEFINE(HAVE_SUNDIALS_SUNCONTEXT, 1, [Define to 1 if SUNDIALS' API is using a SUNContext object.])
if test $octave_have_sundials_compatible_api = no; then
warn_sundials_disabled="SUNDIALS libraries do not provide an API that is compatible with Octave. The solvers ode15i and ode15s will be disabled."
OCTAVE_CONFIGURE_WARNING([warn_sundials_disabled])
@@ -2746,11 +2747,12 @@ AC_DEFUN([OCTAVE_CHECK_SUNDIALS_SUNLINSO
# include <ufsparse/klu.h>
#endif
])
+ ## Check for current KLU function name first.
OCTAVE_CHECK_LIB(sundials_sunlinsolklu, SUNLINSOL_KLU, [],
- [], [SUNKLU], [],
+ [], [SUNLinSol_KLU], [],
[don't use SUNDIALS SUNLINSOL_KLU library, disable ode15i and ode15s sparse Jacobian],
- [AC_CHECK_FUNCS([SUNLinSol_KLU SUNKLU])
- AC_CACHE_CHECK([whether compiling a program that calls SUNKLU works],
+ [AC_CHECK_FUNCS([SUNLinSol_KLU])
+ AC_CACHE_CHECK([whether compiling a program that calls SUNLinSol_KLU works],
[octave_cv_sundials_sunlinsol_klu],
[AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[
#if defined (HAVE_IDA_IDA_H)
@@ -2772,11 +2774,53 @@ AC_DEFUN([OCTAVE_CHECK_SUNDIALS_SUNLINSO
#include <sunlinsol/sunlinsol_klu.h>
#endif
]], [[
- SUNKLU (0, 0);
+ #if defined (HAVE_SUNCONTEXT_CREATE)
+ SUNContext *sunContext;
+ if (SUNContext_Create (NULL, sunContext) < 0)
+ 1/0; // provoke an error
+ SUNLinSol_KLU (0, 0, *sunContext);
+ SUNContext_Free (sunContext);
+ #else
+ SUNLinSol_KLU (0, 0);
+ #endif
]])],
octave_cv_sundials_sunlinsol_klu=yes,
octave_cv_sundials_sunlinsol_klu=no)
])])
+ if test "x$octave_cv_sundials_sunlinsol_klu" = xno; then
+ ## Check for deprecated KLU function name second.
+ OCTAVE_CHECK_LIB(sundials_sunlinsolklu, SUNLINSOL_KLU, [],
+ [], [SUNKLU], [],
+ [don't use SUNDIALS SUNLINSOL_KLU library, disable ode15i and ode15s sparse Jacobian],
+ [AC_CHECK_FUNCS([SUNKLU])
+ AC_CACHE_CHECK([whether compiling a program that calls SUNLinSol_KLU works],
+ [octave_cv_sundials_sunlinsol_klu],
+ [AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[
+ #if defined (HAVE_IDA_IDA_H)
+ #include <ida/ida.h>
+ #endif
+ #if defined (HAVE_KLU_H)
+ #include <klu.h>
+ #endif
+ #if defined (HAVE_KLU_KLU_H)
+ #include <klu/klu.h>
+ #endif
+ #if defined (HAVE_SUITESPARSE_KLU_H)
+ #include <suitesparse/klu.h>
+ #endif
+ #if defined (HAVE_UFPARSE_KLU_H)
+ #include <ufsparse/klu.h>
+ #endif
+ #if defined (HAVE_SUNLINSOL_SUNLINSOL_KLU_H)
+ #include <sunlinsol/sunlinsol_klu.h>
+ #endif
+ ]], [[
+ SUNKLU (0, 0);
+ ]])],
+ octave_cv_sundials_sunlinsol_klu=yes,
+ octave_cv_sundials_sunlinsol_klu=no)
+ ])])
+ fi
if test "x$ac_cv_header_sunlinsol_sunlinsol_klu_h" = xyes \
&& test "x$octave_cv_sundials_sunlinsol_klu" = xyes; then
AC_DEFINE(HAVE_SUNDIALS_SUNLINSOL_KLU, 1,