Merge remote-tracking branch 'origin/master' into windows
[scilab.git] / scilab / modules / differential_equations / src / cpp / differentialequationfunctions.cpp
1 /*
2  *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  *  Copyright (C) 2011 - DIGITEO - Cedric DELAMARRE
4  *
5  * Copyright (C) 2012 - 2016 - Scilab Enterprises
6  *
7  * This file is hereby licensed under the terms of the GNU GPL v2.0,
8  * pursuant to article 5.3.4 of the CeCILL v.2.1.
9  * This file was originally licensed under the terms of the CeCILL v2.1,
10  * and continues to be available under such terms.
11  * For more information, see the COPYING file which you should have received
12  * along with this program.
13  *
14  */
15 /*--------------------------------------------------------------------------*/
16 #include <list>
17
18 #include "string.hxx"
19 #include "double.hxx"
20 #include "differentialequationfunctions.hxx"
21 #include "configvariable.hxx"
22 #include "commentexp.hxx"
23
24 extern "C"
25 {
26 #include "elem_common.h"
27 #include "scifunctions.h"
28 #include "Ex-odedc.h"
29 #include "Ex-ode.h"
30 #include "Ex-daskr.h"
31 #include "localization.h"
32 }
33
34 /*
35 ** differential equation functions
36 ** \{
37 */
38
39 // need the current thread, not the last running thread.
40
41 std::list<DifferentialEquationFunctions*> DifferentialEquation::m_DifferentialEquationFunctions;
42
43 using namespace types;
44 void DifferentialEquation::addDifferentialEquationFunctions(DifferentialEquationFunctions* _deFunction)
45 {
46     m_DifferentialEquationFunctions.push_back(_deFunction);
47 }
48
49 void DifferentialEquation::removeDifferentialEquationFunctions()
50 {
51     m_DifferentialEquationFunctions.pop_back();
52 }
53
54 DifferentialEquationFunctions* DifferentialEquation::getDifferentialEquationFunctions()
55 {
56     return m_DifferentialEquationFunctions.back();
57 }
58
59
60 /*
61 ** \}
62 */
63 /*--------------------------------------------------------------------------*/
64 DifferentialEquationFunctions::DifferentialEquationFunctions(const std::string& callerName)
65 {
66     m_odeYRows      = 0;
67     m_odeYCols      = 0;
68     m_odedcYDSize   = 0;
69     m_odedcFlag     = 0;
70     m_bvodeM        = 0;
71     m_bvodeN        = 0;
72     m_mu            = 0;
73     m_ml            = 0;
74     m_bandedJac     = false;
75
76     m_strCaller = callerName;
77
78     // callable
79     m_pCallFFunction      = NULL;
80     m_pCallJacFunction    = NULL;
81     m_pCallGFunction      = NULL;
82     m_pCallPjacFunction   = NULL;
83     m_pCallPsolFunction   = NULL;
84
85     // function extern
86     m_pStringFFunctionDyn       = NULL;
87     m_pStringJacFunctionDyn     = NULL;
88     m_pStringGFunctionDyn       = NULL;
89     m_pStringPjacFunctionDyn    = NULL;
90     m_pStringPsolFunctionDyn    = NULL;
91
92     // function static
93     m_pStringFFunctionStatic    = NULL;
94     m_pStringJacFunctionStatic  = NULL;
95     m_pStringGFunctionStatic    = NULL;
96     m_pStringPjacFunctionStatic = NULL;
97     m_pStringPsolFunctionStatic = NULL;
98
99     // bvode
100     m_pCallFsubFunction     = NULL;
101     m_pCallDfsubFunction    = NULL;
102     m_pCallGsubFunction     = NULL;
103     m_pCallDgsubFunction    = NULL;
104     m_pCallGuessFunction    = NULL;
105
106     m_pStringFsubFunctionDyn    = NULL;
107     m_pStringDfsubFunctionDyn   = NULL;
108     m_pStringGsubFunctionDyn    = NULL;
109     m_pStringDgsubFunctionDyn   = NULL;
110     m_pStringGuessFunctionDyn   = NULL;
111
112     m_pStringFsubFunctionStatic     = NULL;
113     m_pStringDfsubFunctionStatic    = NULL;
114     m_pStringGsubFunctionStatic     = NULL;
115     m_pStringDgsubFunctionStatic    = NULL;
116     m_pStringGuessFunctionStatic    = NULL;
117
118     // init static functions
119     if (callerName == "ode")
120     {
121         m_staticFunctionMap["arnol"]   = (void*) C2F(arnol);
122         m_staticFunctionMap["fex"]     = (void*) fex;
123         m_staticFunctionMap["fex2"]    = (void*) fex2;
124         m_staticFunctionMap["fex3"]    = (void*) fex3;
125         m_staticFunctionMap["fexab"]   = (void*) fexab;
126         m_staticFunctionMap["loren"]   = (void*) C2F(loren);
127         m_staticFunctionMap["bcomp"]   = (void*) C2F(bcomp);
128         m_staticFunctionMap["lcomp"]   = (void*) C2F(lcomp);
129
130         m_staticFunctionMap["jex"]     = (void*) jex;
131     }
132     else if (callerName == "odedc")
133     {
134         m_staticFunctionMap["fcd"]     = (void*) fcd;
135         m_staticFunctionMap["fcd1"]    = (void*) fcd1;
136         m_staticFunctionMap["fexcd"]   = (void*) fexcd;
137         m_staticFunctionMap["phis"]    = (void*) phis;
138         m_staticFunctionMap["phit"]    = (void*) phit;
139
140         m_staticFunctionMap["jex"]     = (void*) jex;
141     }
142     else if (callerName == "intg")
143     {
144         m_staticFunctionMap["intgex"]  = (void*) C2F(intgex);
145     }
146     else if (callerName == "int2d")
147     {
148         m_staticFunctionMap["int2dex"] = (void*) C2F(int2dex);
149     }
150     else if (callerName == "int3d")
151     {
152         m_staticFunctionMap["int3dex"] = (void*) C2F(int3dex);
153     }
154     else if (callerName == "feval")
155     {
156         m_staticFunctionMap["parab"]   = (void*) C2F(parab);
157         m_staticFunctionMap["parabc"]  = (void*) C2F(parabc);
158     }
159     else if (callerName == "bvode")
160     {
161         m_staticFunctionMap["cndg"]    = (void*) C2F(cndg);
162         m_staticFunctionMap["cng"]     = (void*) C2F(cng);
163         m_staticFunctionMap["cnf"]     = (void*) C2F(cnf);
164         m_staticFunctionMap["cndf"]    = (void*) C2F(cndf);
165         m_staticFunctionMap["cngu"]    = (void*) C2F(cngu);
166     }
167     else if (callerName == "impl")
168     {
169         m_staticFunctionMap["resid"]   = (void*) C2F(resid);  // res
170         m_staticFunctionMap["aplusp"]  = (void*) C2F(aplusp); // adda
171         m_staticFunctionMap["dgbydy"]  = (void*) C2F(dgbydy); // jac
172     }
173     else if (callerName == "dassl" ||
174              callerName == "dasrt" ||
175              callerName == "daskr")
176     {
177         //res
178         m_staticFunctionMap["res1"]    = (void*) C2F(res1);
179         m_staticFunctionMap["res2"]    = (void*) C2F(res2);
180         m_staticFunctionMap["dres1"]   = (void*) C2F(dres1);
181         m_staticFunctionMap["dres2"]   = (void*) C2F(dres2);
182
183         // jac
184         m_staticFunctionMap["jac2"]   = (void*) C2F(jac2);
185         m_staticFunctionMap["djac2"]  = (void*) C2F(djac2);
186         m_staticFunctionMap["djac1"]  = (void*) C2F(djac1);
187
188         //g
189         if (callerName == "dasrt" || callerName == "daskr")
190         {
191             m_staticFunctionMap["gr1"]  = (void*) C2F(gr1);
192             m_staticFunctionMap["gr2"]  = (void*) C2F(gr2);
193         }
194
195         // pjac, psol
196         if (callerName == "daskr")
197         {
198             m_staticFunctionMap["pjac1"]  = (void*) pjac1;
199             m_staticFunctionMap["psol1"]  = (void*) psol1;
200         }
201     }
202 }
203
204 DifferentialEquationFunctions::~DifferentialEquationFunctions()
205 {
206     m_staticFunctionMap.clear();
207 }
208
209 /*------------------------------- public -------------------------------------------*/
210 void DifferentialEquationFunctions::execDasrtG(int* ny, double* t, double* y, int* ng, double* gout, double* rpar, int* ipar)
211 {
212     char errorMsg[256];
213     if (m_pCallGFunction)
214     {
215         callDasrtMacroG(ny, t, y, ng, gout, rpar, ipar);
216     }
217     else if (m_pStringGFunctionDyn)
218     {
219         ConfigVariable::EntryPointStr* func = ConfigVariable::getEntryPoint(m_pStringGFunctionDyn->get(0));
220         if (func == NULL)
221         {
222             sprintf(errorMsg, _("Undefined function '%s'.\n"), m_pStringGFunctionDyn->get(0));
223             throw ast::InternalError(errorMsg);
224         }
225         ((dasrt_g_t)(func->functionPtr))(ny, t, y, ng, gout, rpar, ipar);
226     }
227     else if (m_pStringGFunctionStatic)
228     {
229         ((dasrt_g_t)m_staticFunctionMap[m_pStringGFunctionStatic->get(0)])(ny, t, y, ng, gout, rpar, ipar);
230     }
231     else
232     {
233         sprintf(errorMsg, _("User function '%s' have not been setted.\n"), "g");
234         throw ast::InternalError(errorMsg);
235     }
236 }
237
238 void DifferentialEquationFunctions::execDasslF(double* t, double* y, double* ydot, double* delta, int* ires, double* rpar, int* ipar)
239 {
240     char errorMsg[256];
241     if (m_pCallFFunction)
242     {
243         callDasslMacroF(t, y, ydot, delta, ires, rpar, ipar);
244     }
245     else if (m_pStringFFunctionDyn)
246     {
247         ConfigVariable::EntryPointStr* func = ConfigVariable::getEntryPoint(m_pStringFFunctionDyn->get(0));
248         if (func == NULL)
249         {
250             sprintf(errorMsg, _("Undefined function '%s'.\n"), m_pStringFFunctionDyn->get(0));
251             throw ast::InternalError(errorMsg);
252         }
253         ((dassl_f_t)(func->functionPtr))(t, y, ydot, delta, ires, rpar, ipar);
254     }
255     else if (m_pStringFFunctionStatic)
256     {
257         ((dassl_f_t)m_staticFunctionMap[m_pStringFFunctionStatic->get(0)])(t, y, ydot, delta, ires, rpar, ipar);
258     }
259     else
260     {
261         sprintf(errorMsg, _("User function '%s' have not been setted.\n"), "f");
262         throw ast::InternalError(errorMsg);
263     }
264 }
265
266 void DifferentialEquationFunctions::execDasslJac(double* t, double* y, double* ydot, double* pd, double* cj, double* rpar, int* ipar)
267 {
268     char errorMsg[256];
269     if (m_pCallJacFunction)
270     {
271         callDasslMacroJac(t, y, ydot, pd, cj, rpar, ipar);
272     }
273     else if (m_pStringJacFunctionDyn)
274     {
275         ConfigVariable::EntryPointStr* func = ConfigVariable::getEntryPoint(m_pStringJacFunctionDyn->get(0));
276         if (func == NULL)
277         {
278             sprintf(errorMsg, _("Undefined function '%s'.\n"), m_pStringJacFunctionDyn->get(0));
279             throw ast::InternalError(errorMsg);
280         }
281         ((dassl_jac_t)(func->functionPtr))(t, y, ydot, pd, cj, rpar, ipar);
282     }
283     else if (m_pStringJacFunctionStatic)
284     {
285         ((dassl_jac_t)m_staticFunctionMap[m_pStringJacFunctionStatic->get(0)])(t, y, ydot, pd, cj, rpar, ipar);
286     }
287     else
288     {
289         sprintf(errorMsg, _("User function '%s' have not been setted.\n"), "jacobian");
290         throw ast::InternalError(errorMsg);
291     }
292 }
293
294 void DifferentialEquationFunctions::execDaskrPjac(double* res, int* ires, int* neq, double* t, double* y, double* ydot,
295         double* rewt, double* savr, double* wk, double* h, double* cj,
296         double* wp, int* iwp, int* ier, double* rpar, int* ipar)
297 {
298     char errorMsg[256];
299     if (m_pCallPjacFunction)
300     {
301         callDaskrMacroPjac(res, ires, neq, t, y, ydot, rewt, savr,
302                            wk, h, cj, wp, iwp, ier, rpar, ipar);
303     }
304     else if (m_pStringPjacFunctionDyn)
305     {
306         ConfigVariable::EntryPointStr* func = ConfigVariable::getEntryPoint(m_pStringPjacFunctionDyn->get(0));
307         if (func == NULL)
308         {
309             sprintf(errorMsg, _("Undefined function '%s'.\n"), m_pStringPjacFunctionDyn->get(0));
310             throw ast::InternalError(errorMsg);
311         }
312         ((daskr_pjac_t)(func->functionPtr))(res, ires, neq, t, y, ydot, rewt, savr,
313                                             wk, h, cj, wp, iwp, ier, rpar, ipar);
314     }
315     else if (m_pStringPjacFunctionStatic)
316     {
317         ((daskr_pjac_t)m_staticFunctionMap[m_pStringPjacFunctionStatic->get(0)])(res, ires, neq, t, y, ydot, rewt, savr,
318                 wk, h, cj, wp, iwp, ier, rpar, ipar);
319     }
320     else
321     {
322         sprintf(errorMsg, _("User function '%s' have not been setted.\n"), "pjac");
323         throw ast::InternalError(errorMsg);
324     }
325 }
326
327 void DifferentialEquationFunctions::execDaskrPsol(int* neq, double* t, double* y, double* ydot, double* savr, double* wk,
328         double* cj, double* wght, double* wp, int* iwp, double* b, double* eplin,
329         int* ier, double* rpar, int* ipar)
330 {
331     char errorMsg[256];
332     if (m_pCallPsolFunction)
333     {
334         callDaskrMacroPsol(neq, t, y, ydot, savr, wk, cj, wght,
335                            wp, iwp, b, eplin, ier, rpar, ipar);
336     }
337     else if (m_pStringPsolFunctionDyn)
338     {
339         ConfigVariable::EntryPointStr* func = ConfigVariable::getEntryPoint(m_pStringPsolFunctionDyn->get(0));
340         if (func == NULL)
341         {
342             sprintf(errorMsg, _("Undefined function '%s'.\n"), m_pStringPsolFunctionDyn->get(0));
343             throw ast::InternalError(errorMsg);
344         }
345         ((daskr_psol_t)(func->functionPtr))(neq, t, y, ydot, savr, wk, cj, wght,
346                                             wp, iwp, b, eplin, ier, rpar, ipar);
347     }
348     else if (m_pStringPsolFunctionStatic)
349     {
350         ((daskr_psol_t)m_staticFunctionMap[m_pStringPsolFunctionStatic->get(0)])(neq, t, y, ydot, savr, wk, cj, wght,
351                 wp, iwp, b, eplin, ier, rpar, ipar);
352     }
353     else
354     {
355         sprintf(errorMsg, _("User function '%s' have not been setted.\n"), "psol");
356         throw ast::InternalError(errorMsg);
357     }
358 }
359
360 void DifferentialEquationFunctions::execImplF(int* neq, double* t, double* y, double* s, double* r, int* ires)
361 {
362     char errorMsg[256];
363     if (m_pCallFFunction)
364     {
365         callImplMacroF(neq, t, y, s, r, ires);
366     }
367     else if (m_pStringFFunctionDyn)
368     {
369         ConfigVariable::EntryPointStr* func = ConfigVariable::getEntryPoint(m_pStringFFunctionDyn->get(0));
370         if (func == NULL)
371         {
372             sprintf(errorMsg, _("Undefined function '%s'.\n"), m_pStringFFunctionDyn->get(0));
373             throw ast::InternalError(errorMsg);
374         }
375         ((impl_f_t)(func->functionPtr))(neq, t, y, s, r, ires);
376     }
377     else if (m_pStringFFunctionStatic)
378     {
379         ((impl_f_t)m_staticFunctionMap[m_pStringFFunctionStatic->get(0)])(neq, t, y, s, r, ires);
380     }
381     else
382     {
383         sprintf(errorMsg, _("User function '%s' have not been setted.\n"), "f");
384         throw ast::InternalError(errorMsg);
385     }
386 }
387
388 void DifferentialEquationFunctions::execImplG(int* neq, double* t, double* y, double* ml, double* mu, double* p, int* nrowp)
389 {
390     char errorMsg[256];
391     if (m_pCallGFunction)
392     {
393         callImplMacroG(neq, t, y, ml, mu, p, nrowp);
394     }
395     else if (m_pStringGFunctionDyn)
396     {
397         ConfigVariable::EntryPointStr* func = ConfigVariable::getEntryPoint(m_pStringGFunctionDyn->get(0));
398         if (func == NULL)
399         {
400             sprintf(errorMsg, _("Undefined function '%s'.\n"), m_pStringGFunctionDyn->get(0));
401             throw ast::InternalError(errorMsg);
402         }
403         ((impl_g_t)(func->functionPtr))(neq, t, y, ml, mu, p, nrowp);
404     }
405     else if (m_pStringGFunctionStatic)
406     {
407         ((impl_g_t)m_staticFunctionMap[m_pStringGFunctionStatic->get(0)])(neq, t, y, ml, mu, p, nrowp);
408     }
409     else
410     {
411         sprintf(errorMsg, _("User function '%s' have not been setted.\n"), "g");
412         throw ast::InternalError(errorMsg);
413     }
414 }
415
416 void DifferentialEquationFunctions::execImplJac(int* neq, double* t, double* y, double* s, double* ml, double* mu, double* p, int* nrowp)
417 {
418     char errorMsg[256];
419     if (m_pCallJacFunction)
420     {
421         callImplMacroJac(neq, t, y, s, ml, mu, p, nrowp);
422     }
423     else if (m_pStringJacFunctionDyn)
424     {
425         ConfigVariable::EntryPointStr* func = ConfigVariable::getEntryPoint(m_pStringJacFunctionDyn->get(0));
426         if (func == NULL)
427         {
428             sprintf(errorMsg, _("Undefined function '%s'.\n"), m_pStringJacFunctionDyn->get(0));
429             throw ast::InternalError(errorMsg);
430         }
431         ((impl_jac_t)(func->functionPtr))(neq, t, y, s, ml, mu, p, nrowp);
432     }
433     else if (m_pStringJacFunctionStatic)
434     {
435         ((impl_jac_t)m_staticFunctionMap[m_pStringJacFunctionStatic->get(0)])(neq, t, y, s, ml, mu, p, nrowp);
436     }
437     else
438     {
439         sprintf(errorMsg, _("User function '%s' have not been setted.\n"), "jacobian");
440         throw ast::InternalError(errorMsg);
441     }
442 }
443
444 void DifferentialEquationFunctions::execBvodeGuess(double *x, double *z, double *d)
445 {
446     char errorMsg[256];
447     if (m_pCallGuessFunction)
448     {
449         callBvodeMacroGuess(x, z, d);
450     }
451     else if (m_pStringGuessFunctionDyn)
452     {
453         ConfigVariable::EntryPointStr* func = ConfigVariable::getEntryPoint(m_pStringGuessFunctionDyn->get(0));
454         if (func == NULL)
455         {
456             sprintf(errorMsg, _("Undefined function '%s'.\n"), m_pStringGuessFunctionDyn->get(0));
457             throw ast::InternalError(errorMsg);
458         }
459         ((bvode_ddd_t)(func->functionPtr))(x, z, d);
460     }
461     else if (m_pStringGuessFunctionStatic)
462     {
463         ((bvode_ddd_t)m_staticFunctionMap[m_pStringGuessFunctionStatic->get(0)])(x, z, d);
464     }
465     else
466     {
467         sprintf(errorMsg, _("User function '%s' have not been setted.\n"), "guess");
468         throw ast::InternalError(errorMsg);
469     }
470 }
471
472 void DifferentialEquationFunctions::execBvodeDfsub(double *x, double *z, double *d)
473 {
474     char errorMsg[256];
475     if (m_pCallDfsubFunction)
476     {
477         callBvodeMacroDfsub(x, z, d);
478     }
479     else if (m_pStringDfsubFunctionDyn)
480     {
481         ConfigVariable::EntryPointStr* func = ConfigVariable::getEntryPoint(m_pStringDfsubFunctionDyn->get(0));
482         if (func == NULL)
483         {
484             sprintf(errorMsg, _("Undefined function '%s'.\n"), m_pStringDfsubFunctionDyn->get(0));
485             throw ast::InternalError(errorMsg);
486         }
487         ((bvode_ddd_t)(func->functionPtr))(x, z, d);
488     }
489     else if (m_pStringDfsubFunctionStatic)// function static
490     {
491         ((bvode_ddd_t)m_staticFunctionMap[m_pStringDfsubFunctionStatic->get(0)])(x, z, d);
492     }
493     else
494     {
495         sprintf(errorMsg, _("User function '%s' have not been setted.\n"), "fsub");
496         throw ast::InternalError(errorMsg);
497     }
498 }
499
500 void DifferentialEquationFunctions::execBvodeFsub(double *x, double *z, double *d)
501 {
502     char errorMsg[256];
503     if (m_pCallFsubFunction)
504     {
505         callBvodeMacroFsub(x, z, d);
506     }
507     else if (m_pStringFsubFunctionDyn)
508     {
509         ConfigVariable::EntryPointStr* func = ConfigVariable::getEntryPoint(m_pStringFsubFunctionDyn->get(0));
510         if (func == NULL)
511         {
512             sprintf(errorMsg, _("Undefined function '%s'.\n"), m_pStringFsubFunctionDyn->get(0));
513             throw ast::InternalError(errorMsg);
514         }
515         ((bvode_ddd_t)(func->functionPtr))(x, z, d);
516     }
517     else if (m_pStringFsubFunctionStatic) // function static
518     {
519         ((bvode_ddd_t)m_staticFunctionMap[m_pStringFsubFunctionStatic->get(0)])(x, z, d);
520     }
521     else
522     {
523         sprintf(errorMsg, _("User function '%s' have not been setted.\n"), "fsub");
524         throw ast::InternalError(errorMsg);
525     }
526 }
527
528 void DifferentialEquationFunctions::execBvodeDgsub(int *i, double *z, double *g)
529 {
530     char errorMsg[256];
531     if (m_pCallDgsubFunction)
532     {
533         callBvodeMacroDgsub(i, z, g);
534     }
535     else if (m_pStringDgsubFunctionDyn)
536     {
537         ConfigVariable::EntryPointStr* func = ConfigVariable::getEntryPoint(m_pStringDgsubFunctionDyn->get(0));
538         if (func == NULL)
539         {
540             sprintf(errorMsg, _("Undefined function '%s'.\n"), m_pStringDgsubFunctionDyn->get(0));
541             throw ast::InternalError(errorMsg);
542         }
543         ((bvode_idd_t)(func->functionPtr))(i, z, g);
544     }
545     else if (m_pStringDgsubFunctionStatic) // function static
546     {
547         ((bvode_idd_t)m_staticFunctionMap[m_pStringDgsubFunctionStatic->get(0)])(i, z, g);
548     }
549     else
550     {
551         sprintf(errorMsg, _("User function '%s' have not been setted.\n"), "gsub");
552         throw ast::InternalError(errorMsg);
553     }
554 }
555
556 void DifferentialEquationFunctions::execBvodeGsub(int *i, double *z, double *g)
557 {
558     char errorMsg[256];
559     if (m_pCallGsubFunction)
560     {
561         callBvodeMacroGsub(i, z, g);
562     }
563     else if (m_pStringGsubFunctionDyn)
564     {
565         ConfigVariable::EntryPointStr* func = ConfigVariable::getEntryPoint(m_pStringGsubFunctionDyn->get(0));
566         if (func == NULL)
567         {
568             sprintf(errorMsg, _("Undefined function '%s'.\n"), m_pStringGsubFunctionDyn->get(0));
569             throw ast::InternalError(errorMsg);
570         }
571         ((bvode_idd_t)(func->functionPtr))(i, z, g);
572     }
573     else if (m_pStringGsubFunctionStatic) // function static
574     {
575         ((bvode_idd_t)m_staticFunctionMap[m_pStringGsubFunctionStatic->get(0)])(i, z, g);
576     }
577     else
578     {
579         sprintf(errorMsg, _("User function '%s' have not been setted.\n"), "gsub");
580         throw ast::InternalError(errorMsg);
581     }
582 }
583
584 void DifferentialEquationFunctions::execFevalF(int *nn, double *x1, double *x2, double *xres, int *itype)
585 {
586     char errorMsg[256];
587     if (m_pCallFFunction)
588     {
589         callFevalMacroF(nn, x1, x2, xres, itype);
590     }
591     else if (m_pStringFFunctionDyn)
592     {
593         ConfigVariable::EntryPointStr* func = ConfigVariable::getEntryPoint(m_pStringFFunctionDyn->get(0));
594         if (func == NULL)
595         {
596             sprintf(errorMsg, _("Undefined function '%s'.\n"), m_pStringFFunctionDyn->get(0));
597             throw ast::InternalError(errorMsg);
598         }
599
600         ((feval_f_t)(func->functionPtr))(nn, x1, x2, xres, itype);
601     }
602     else if (m_pStringFFunctionStatic) // function static
603     {
604         ((feval_f_t)m_staticFunctionMap[m_pStringFFunctionStatic->get(0)])(nn, x1, x2, xres, itype);
605     }
606     else
607     {
608         sprintf(errorMsg, _("User function '%s' have not been setted.\n"), "f");
609         throw ast::InternalError(errorMsg);
610     }
611 }
612
613 void DifferentialEquationFunctions::execInt3dF(double* x, int* numfun, double* funvls)
614 {
615     char errorMsg[256];
616     if (m_pCallFFunction)
617     {
618         callInt3dMacroF(x, numfun, funvls);
619     }
620     else if (m_pStringFFunctionDyn)
621     {
622         ConfigVariable::EntryPointStr* func = ConfigVariable::getEntryPoint(m_pStringFFunctionDyn->get(0));
623         if (func == NULL)
624         {
625             sprintf(errorMsg, _("Undefined function '%s'.\n"), m_pStringFFunctionDyn->get(0));
626             throw ast::InternalError(errorMsg);
627         }
628         ((int3d_f_t)(func->functionPtr))(x, numfun, funvls);
629     }
630     else if (m_pStringFFunctionStatic) // function static
631     {
632         ((int3d_f_t)m_staticFunctionMap[m_pStringFFunctionStatic->get(0)])(x, numfun, funvls);
633     }
634     else
635     {
636         sprintf(errorMsg, _("User function '%s' have not been setted.\n"), "f");
637         throw ast::InternalError(errorMsg);
638     }
639 }
640
641 double DifferentialEquationFunctions::execInt2dF(double* x, double* y)
642 {
643     char errorMsg[256];
644     if (m_pCallFFunction)
645     {
646         return callInt2dMacroF(x, y);
647     }
648     else if (m_pStringFFunctionDyn)
649     {
650         ConfigVariable::EntryPointStr* func = ConfigVariable::getEntryPoint(m_pStringFFunctionDyn->get(0));
651         if (func == NULL)
652         {
653             sprintf(errorMsg, _("Undefined function '%s'.\n"), m_pStringFFunctionDyn->get(0));
654             throw ast::InternalError(errorMsg);
655         }
656         return ((int2d_f_t)(func->functionPtr))(x, y);
657     }
658     else if (m_pStringFFunctionStatic) // function static
659     {
660         return ((int2d_f_t)m_staticFunctionMap[m_pStringFFunctionStatic->get(0)])(x, y);
661     }
662     else
663     {
664         sprintf(errorMsg, _("User function '%s' have not been setted.\n"), "f");
665         throw ast::InternalError(errorMsg);
666     }
667 }
668
669 double DifferentialEquationFunctions::execIntgF(double* x)
670 {
671     char errorMsg[256];
672     if (m_pCallFFunction)
673     {
674         return callIntgMacroF(x);
675     }
676     else if (m_pStringFFunctionDyn)
677     {
678         ConfigVariable::EntryPointStr* func = ConfigVariable::getEntryPoint(m_pStringFFunctionDyn->get(0));
679         if (func == NULL)
680         {
681             sprintf(errorMsg, _("Undefined function '%s'.\n"), m_pStringFFunctionDyn->get(0));
682             throw ast::InternalError(errorMsg);
683         }
684         return ((intg_f_t)(func->functionPtr))(x);
685     }
686     else if (m_pStringFFunctionStatic) // function static
687     {
688         return ((intg_f_t)m_staticFunctionMap[m_pStringFFunctionStatic->get(0)])(x);
689     }
690     else
691     {
692         sprintf(errorMsg, _("User function '%s' have not been setted.\n"), "f");
693         throw ast::InternalError(errorMsg);
694     }
695 }
696
697 void DifferentialEquationFunctions::execOdeF(int* n, double* t, double* y, double* yout)
698 {
699     char errorMsg[256];
700     if (m_pCallFFunction)
701     {
702         callOdeMacroF(n, t, y, yout);
703     }
704     else if (m_pStringFFunctionDyn)
705     {
706         ConfigVariable::EntryPointStr* func = ConfigVariable::getEntryPoint(m_pStringFFunctionDyn->get(0));
707         if (func == NULL)
708         {
709             sprintf(errorMsg, _("Undefined function '%s'.\n"), m_pStringFFunctionDyn->get(0));
710             throw ast::InternalError(errorMsg);
711         }
712
713         if (m_strCaller == "ode")
714         {
715             ((ode_f_t)(func->functionPtr))(n, t, y, yout);
716         }
717         else
718         {
719             ((odedc_f_t)(func->functionPtr))(&m_odedcFlag, n, &m_odedcYDSize, t, y, yout);
720         }
721     }
722     else if (m_pStringFFunctionStatic) // function static
723     {
724         if (m_strCaller == "ode")
725         {
726             ((ode_f_t)m_staticFunctionMap[m_pStringFFunctionStatic->get(0)])(n, t, y, yout);
727         }
728         else // if (m_strCaller == "odedc")
729         {
730             ((odedc_f_t)m_staticFunctionMap[m_pStringFFunctionStatic->get(0)])(&m_odedcFlag, n, &m_odedcYDSize, t, y, yout);
731         }
732     }
733     else
734     {
735         sprintf(errorMsg, _("User function '%s' have not been setted.\n"), "f");
736         throw ast::InternalError(errorMsg);
737     }
738 }
739
740 void DifferentialEquationFunctions::execFunctionJac(int *n, double *t, double *y, int *ml, int *mu, double *J, int *nrpd)
741 {
742     char errorMsg[256];
743     if (m_pCallJacFunction)
744     {
745         callMacroJac(n, t, y, ml, mu, J, nrpd);
746     }
747     else if (m_pStringJacFunctionDyn)
748     {
749         ConfigVariable::EntryPointStr* func = ConfigVariable::getEntryPoint(m_pStringJacFunctionDyn->get(0));
750         if (func == NULL)
751         {
752             sprintf(errorMsg, _("Undefined function '%s'.\n"), m_pStringJacFunctionDyn->get(0));
753             throw ast::InternalError(errorMsg);
754         }
755         ((func_jac_t)(func->functionPtr))(n, t, y, ml, mu, J, nrpd);
756     }
757     else if (m_pStringJacFunctionStatic) // function static
758     {
759         ((func_jac_t)m_staticFunctionMap[m_pStringJacFunctionStatic->get(0)])(n, t, y, ml, mu, J, nrpd);
760     }
761     else
762     {
763         sprintf(errorMsg, _("User function '%s' have not been setted.\n"), "jacobian");
764         throw ast::InternalError(errorMsg);
765     }
766 }
767
768 void DifferentialEquationFunctions::execFunctionG(int* n, double* t, double* y, int* ng, double* gout)
769 {
770     char errorMsg[256];
771     if (m_pCallGFunction)
772     {
773         callMacroG(n, t, y, ng, gout);
774     }
775     else if (m_pStringGFunctionDyn)
776     {
777         ConfigVariable::EntryPointStr* func = ConfigVariable::getEntryPoint(m_pStringGFunctionDyn->get(0));
778         if (func == NULL)
779         {
780             sprintf(errorMsg, _("Undefined function '%s'.\n"), m_pStringGFunctionDyn->get(0));
781             throw ast::InternalError(errorMsg);
782         }
783         ((func_g_t)(func->functionPtr))(n, t, y, ng, gout);
784     }
785     else if (m_pStringGFunctionStatic)// function static
786     {
787         ((func_g_t)m_staticFunctionMap[m_pStringGFunctionStatic->get(0)])(n, t, y, ng, gout);
788     }
789     else
790     {
791         sprintf(errorMsg, _("User function '%s' have not been setted.\n"), "g");
792         throw ast::InternalError(errorMsg);
793     }
794 }
795
796 //*** setter ***
797 // set rows cols
798 void DifferentialEquationFunctions::setOdeYRows(int rows)
799 {
800     m_odeYRows = rows;
801 }
802
803 void DifferentialEquationFunctions::setOdeYCols(int cols)
804 {
805     m_odeYCols = cols;
806 }
807
808 // set odedc yd size
809 void DifferentialEquationFunctions::setOdedcYDSize(int size)
810 {
811     m_odedcYDSize = size;
812 }
813
814 // set odedc flag
815 void DifferentialEquationFunctions::setOdedcFlag()
816 {
817     m_odedcFlag = 1;
818 }
819
820 // reset odedc flag
821 void DifferentialEquationFunctions::resetOdedcFlag()
822 {
823     m_odedcFlag = 0;
824 }
825
826 void DifferentialEquationFunctions::setBvodeM(int _m)
827 {
828     m_bvodeM = _m;
829 }
830
831 void DifferentialEquationFunctions::setBvodeN(int _n)
832 {
833     m_bvodeN = _n;
834 }
835
836 //set function f, jac, g, psol, pjac as types::Callable
837 void DifferentialEquationFunctions::setFFunction(types::Callable* _odeFFunc)
838 {
839     m_pCallFFunction = _odeFFunc;
840 }
841
842 void DifferentialEquationFunctions::setJacFunction(types::Callable* _odeJacFunc)
843 {
844     m_pCallJacFunction = _odeJacFunc;
845 }
846
847 void DifferentialEquationFunctions::setGFunction(types::Callable* _odeGFunc)
848 {
849     m_pCallGFunction = _odeGFunc;
850 }
851
852 void DifferentialEquationFunctions::setPsolFunction(types::Callable* _pSolFunc)
853 {
854     m_pCallPsolFunction = _pSolFunc;
855 }
856
857 void DifferentialEquationFunctions::setPjacFunction(types::Callable* _pJacFunc)
858 {
859     m_pCallPjacFunction = _pJacFunc;
860 }
861
862 //set function f, jac, g, psol, pjac as types::String
863 bool DifferentialEquationFunctions::setFFunction(types::String* _odeFFunc)
864 {
865     if (ConfigVariable::getEntryPoint(_odeFFunc->get(0)))
866     {
867         m_pStringFFunctionDyn = _odeFFunc;
868         return true;
869     }
870     else
871     {
872         if (m_staticFunctionMap.find(_odeFFunc->get(0)) != m_staticFunctionMap.end())
873         {
874             m_pStringFFunctionStatic = _odeFFunc;
875             return true;
876         }
877         return false;
878     }
879 }
880
881 bool DifferentialEquationFunctions::setJacFunction(types::String* _odeJacFunc)
882 {
883     if (ConfigVariable::getEntryPoint(_odeJacFunc->get(0)))
884     {
885         m_pStringJacFunctionDyn = _odeJacFunc;
886         return true;
887     }
888     else
889     {
890         if (m_staticFunctionMap.find(_odeJacFunc->get(0)) != m_staticFunctionMap.end())
891         {
892             m_pStringJacFunctionStatic = _odeJacFunc;
893             return true;
894         }
895         return false;
896     }
897 }
898
899 bool DifferentialEquationFunctions::setGFunction(types::String* _odeGFunc)
900 {
901     if (ConfigVariable::getEntryPoint(_odeGFunc->get(0)))
902     {
903         m_pStringGFunctionDyn = _odeGFunc;
904         return true;
905     }
906     else
907     {
908         if (m_staticFunctionMap.find(_odeGFunc->get(0)) != m_staticFunctionMap.end())
909         {
910             m_pStringGFunctionStatic = _odeGFunc;
911             return true;
912         }
913         return false;
914     }
915 }
916
917 bool DifferentialEquationFunctions::setPsolFunction(types::String* _pSolFunc)
918 {
919     if (ConfigVariable::getEntryPoint(_pSolFunc->get(0)))
920     {
921         m_pStringPsolFunctionDyn = _pSolFunc;
922         return true;
923     }
924     else
925     {
926         if (m_staticFunctionMap.find(_pSolFunc->get(0)) != m_staticFunctionMap.end())
927         {
928             m_pStringPsolFunctionStatic = _pSolFunc;
929             return true;
930         }
931         return false;
932     }
933 }
934
935 bool DifferentialEquationFunctions::setPjacFunction(types::String* _pJacFunc)
936 {
937     if (ConfigVariable::getEntryPoint(_pJacFunc->get(0)))
938     {
939         m_pStringPjacFunctionDyn = _pJacFunc;
940         return true;
941     }
942     else
943     {
944         if (m_staticFunctionMap.find(_pJacFunc->get(0)) != m_staticFunctionMap.end())
945         {
946             m_pStringPjacFunctionStatic = _pJacFunc;
947             return true;
948         }
949         return false;
950     }
951 }
952
953 // set args for f, jac, g, pjac and psol functions
954 void DifferentialEquationFunctions::setFArgs(types::InternalType* _odeFArg)
955 {
956     m_FArgs.push_back(_odeFArg);
957 }
958
959 void DifferentialEquationFunctions::setJacArgs(types::InternalType* _odeJacArg)
960 {
961     m_JacArgs.push_back(_odeJacArg);
962 }
963
964 void DifferentialEquationFunctions::setGArgs(types::InternalType* _odeGArg)
965 {
966     m_odeGArgs.push_back(_odeGArg);
967 }
968
969 void DifferentialEquationFunctions::setPsolArgs(types::InternalType* _pSolArg)
970 {
971     m_pSolArgs.push_back(_pSolArg);
972 }
973
974 void DifferentialEquationFunctions::setPjacArgs(types::InternalType* _pJacArg)
975 {
976     m_pJacArgs.push_back(_pJacArg);
977 }
978
979 // bvode set function as types::Callable gsub, dgsub, fsub, dfsub, guess
980 void DifferentialEquationFunctions::setGsubFunction(types::Callable* _func)
981 {
982     m_pCallGsubFunction = _func;
983 }
984
985 void DifferentialEquationFunctions::setDgsubFunction(types::Callable* _func)
986 {
987     m_pCallDgsubFunction = _func;
988 }
989
990 void DifferentialEquationFunctions::setFsubFunction(types::Callable* _func)
991 {
992     m_pCallFsubFunction = _func;
993 }
994
995 void DifferentialEquationFunctions::setDfsubFunction(types::Callable* _func)
996 {
997     m_pCallDfsubFunction = _func;
998 }
999
1000 void DifferentialEquationFunctions::setGuessFunction(types::Callable* _func)
1001 {
1002     m_pCallGuessFunction = _func;
1003 }
1004
1005 // bvode set function as types::String gsub, dgsub, fsub, dfsub, guess
1006 bool DifferentialEquationFunctions::setGsubFunction(types::String* _func)
1007 {
1008     if (ConfigVariable::getEntryPoint(_func->get(0)))
1009     {
1010         m_pStringGsubFunctionDyn = _func;
1011         return true;
1012     }
1013     else
1014     {
1015         if (m_staticFunctionMap.find(_func->get(0)) != m_staticFunctionMap.end())
1016         {
1017             m_pStringGsubFunctionStatic = _func;
1018             return true;
1019         }
1020         return false;
1021     }
1022 }
1023
1024 bool DifferentialEquationFunctions::setDgsubFunction(types::String* _func)
1025 {
1026     if (ConfigVariable::getEntryPoint(_func->get(0)))
1027     {
1028         m_pStringDgsubFunctionDyn = _func;
1029         return true;
1030     }
1031     else
1032     {
1033         if (m_staticFunctionMap.find(_func->get(0)) != m_staticFunctionMap.end())
1034         {
1035             m_pStringDgsubFunctionStatic = _func;
1036             return true;
1037         }
1038         return false;
1039     }
1040 }
1041
1042 bool DifferentialEquationFunctions::setFsubFunction(types::String* _func)
1043 {
1044     if (ConfigVariable::getEntryPoint(_func->get(0)))
1045     {
1046         m_pStringFsubFunctionDyn = _func;
1047         return true;
1048     }
1049     else
1050     {
1051         if (m_staticFunctionMap.find(_func->get(0)) != m_staticFunctionMap.end())
1052         {
1053             m_pStringFsubFunctionStatic = _func;
1054             return true;
1055         }
1056         return false;
1057     }
1058 }
1059
1060 bool DifferentialEquationFunctions::setDfsubFunction(types::String* _func)
1061 {
1062     if (ConfigVariable::getEntryPoint(_func->get(0)))
1063     {
1064         m_pStringDfsubFunctionDyn = _func;
1065         return true;
1066     }
1067     else
1068     {
1069         if (m_staticFunctionMap.find(_func->get(0)) != m_staticFunctionMap.end())
1070         {
1071             m_pStringDfsubFunctionStatic = _func;
1072             return true;
1073         }
1074         return false;
1075     }
1076 }
1077
1078 bool DifferentialEquationFunctions::setGuessFunction(types::String* _func)
1079 {
1080     if (ConfigVariable::getEntryPoint(_func->get(0)))
1081     {
1082         m_pStringGuessFunctionDyn = _func;
1083         return true;
1084     }
1085     else
1086     {
1087         if (m_staticFunctionMap.find(_func->get(0)) != m_staticFunctionMap.end())
1088         {
1089             m_pStringGuessFunctionStatic = _func;
1090             return true;
1091         }
1092         return false;
1093     }
1094 }
1095
1096 // bvode set set args for gsub, dgsub, fsub, dfsub, guess functions
1097 void DifferentialEquationFunctions::setGsubArgs(types::InternalType* _arg)
1098 {
1099     m_GsubArgs.push_back(_arg);
1100 }
1101
1102 void DifferentialEquationFunctions::setDgsubArgs(types::InternalType* _arg)
1103 {
1104     m_DgsubArgs.push_back(_arg);
1105 }
1106
1107 void DifferentialEquationFunctions::setFsubArgs(types::InternalType* _arg)
1108 {
1109     m_FsubArgs.push_back(_arg);
1110 }
1111
1112 void DifferentialEquationFunctions::setDfsubArgs(types::InternalType* _arg)
1113 {
1114     m_DfsubArgs.push_back(_arg);
1115 }
1116
1117 void DifferentialEquationFunctions::setGuessArgs(types::InternalType* _arg)
1118 {
1119     m_GuessArgs.push_back(_arg);
1120 }
1121
1122 // set mu and ml
1123 void DifferentialEquationFunctions::setMu(int mu)
1124 {
1125     m_bandedJac = true;
1126     m_mu = mu;
1127 }
1128
1129 void DifferentialEquationFunctions::setMl(int ml)
1130 {
1131     m_bandedJac = true;
1132     m_ml = ml;
1133 }
1134
1135 //*** getter ***
1136 // get y rows cols
1137 int DifferentialEquationFunctions::getOdeYRows()
1138 {
1139     return m_odeYRows;
1140 }
1141
1142 int DifferentialEquationFunctions::getOdeYCols()
1143 {
1144     return m_odeYCols;
1145 }
1146
1147 // get odedc yd size
1148 int DifferentialEquationFunctions::getOdedcYDSize()
1149 {
1150     return m_odedcYDSize;
1151 }
1152
1153 // get odedc flag
1154 int DifferentialEquationFunctions::getOdedcFlag()
1155 {
1156     return m_odedcFlag;
1157 }
1158
1159 /*------------------------------- private -------------------------------------------*/
1160 void DifferentialEquationFunctions::callOdeMacroF(int* n, double* t, double* y, double* ydot)
1161 {
1162     char errorMsg[256];
1163     int one         = 1;
1164     int iRetCount   = 1;
1165
1166     types::typed_list in;
1167     types::typed_list out;
1168     types::optional_list opt;
1169
1170     types::Double* pDblY    = NULL;
1171     types::Double* pDblYC   = NULL;
1172     types::Double* pDblYD   = NULL;
1173     types::Double* pDblFlag = NULL;
1174
1175     // create input args
1176     types::Double* pDblT = new types::Double(*t);
1177     pDblT->IncreaseRef();
1178
1179     if (m_odedcYDSize) // odedc
1180     {
1181         pDblYC = new types::Double(*n, 1);
1182         pDblYC->set(y);
1183         pDblYC->IncreaseRef();
1184         pDblYD = new types::Double(m_odedcYDSize, 1);
1185         pDblYD->set(y + *n);
1186         pDblYD->IncreaseRef();
1187         pDblFlag = new types::Double(m_odedcFlag);
1188         pDblFlag->IncreaseRef();
1189     }
1190     else // ode
1191     {
1192         pDblY = new types::Double(m_odeYRows, m_odeYCols);
1193         pDblY->set(y);
1194         pDblY->IncreaseRef();
1195     }
1196
1197     // push_back
1198     in.push_back(pDblT);
1199     if (m_odedcYDSize) // odedc
1200     {
1201         in.push_back(pDblYC);
1202         in.push_back(pDblYD);
1203         in.push_back(pDblFlag);
1204     }
1205     else
1206     {
1207         in.push_back(pDblY);
1208     }
1209
1210     for (int i = 0; i < (int)m_FArgs.size(); i++)
1211     {
1212         m_FArgs[i]->IncreaseRef();
1213         in.push_back(m_FArgs[i]);
1214     }
1215
1216     try
1217     {
1218         // new std::string("") is delete in destructor of ast::CommentExp
1219         m_pCallFFunction->invoke(in, opt, iRetCount, out, ast::CommentExp(Location(), new std::string("")));
1220     }
1221     catch (const ast::InternalError& ie)
1222     {
1223         for (int i = 0; i < (int)m_FArgs.size(); i++)
1224         {
1225             m_FArgs[i]->DecreaseRef();
1226         }
1227
1228         throw ie;
1229     }
1230
1231     for (int i = 0; i < (int)m_FArgs.size(); i++)
1232     {
1233         m_FArgs[i]->DecreaseRef();
1234     }
1235
1236     if (out.size() != iRetCount)
1237     {
1238         const char* pstrName = m_pCallFFunction->getName().c_str();
1239         sprintf(errorMsg, _("%s: Wrong number of output argument(s): %d expected.\n"), pstrName, iRetCount);
1240         throw ast::InternalError(errorMsg);
1241     }
1242
1243     out[0]->IncreaseRef();
1244
1245     pDblT->DecreaseRef();
1246     if (pDblT->isDeletable())
1247     {
1248         delete pDblT;
1249     }
1250
1251     if (m_odedcYDSize) // odedc
1252     {
1253         pDblYC->DecreaseRef();
1254         if (pDblYC->isDeletable())
1255         {
1256             delete pDblYC;
1257         }
1258         pDblYD->DecreaseRef();
1259         if (pDblYD->isDeletable())
1260         {
1261             delete pDblYD;
1262         }
1263         pDblFlag->DecreaseRef();
1264         if (pDblFlag->isDeletable())
1265         {
1266             delete pDblFlag;
1267         }
1268     }
1269     else
1270     {
1271         pDblY->DecreaseRef();
1272         if (pDblY->isDeletable())
1273         {
1274             delete pDblY;
1275         }
1276     }
1277
1278     out[0]->DecreaseRef();
1279
1280     if (out[0]->isDouble() == false)
1281     {
1282         const char* pstrName = m_pCallFFunction->getName().c_str();
1283         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Real matrix expected.\n"), pstrName, 1);
1284         throw ast::InternalError(errorMsg);
1285     }
1286     types::Double* pDblOut = out[0]->getAs<types::Double>();
1287     if (pDblOut->isComplex())
1288     {
1289         const char* pstrName = m_pCallFFunction->getName().c_str();
1290         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Real matrix expected.\n"), pstrName, 1);
1291         throw ast::InternalError(errorMsg);
1292     }
1293
1294     if (m_odedcFlag && m_odedcYDSize)
1295     {
1296         C2F(dcopy)(&m_odedcYDSize, pDblOut->get(), &one, ydot, &one);
1297     }
1298     else
1299     {
1300         C2F(dcopy)(n, pDblOut->get(), &one, ydot, &one);
1301     }
1302
1303     if (out[0]->isDeletable())
1304     {
1305         delete out[0];
1306     }
1307 }
1308
1309 void DifferentialEquationFunctions::callMacroJac(int* n, double* t, double* y, int* ml, int* mu, double* J, int* nrpd)
1310 {
1311     char errorMsg[256];
1312     int iRetCount   = 1;
1313     int one         = 1;
1314     int iMaxSize    = (*n) * (*nrpd);
1315
1316     types::typed_list in;
1317     types::typed_list out;
1318     types::optional_list opt;
1319
1320     types::Double* pDblY = new types::Double(m_odeYRows, m_odeYCols);
1321     pDblY->set(y);
1322     types::Double* pDblT = new types::Double(*t);
1323
1324     pDblT->IncreaseRef();
1325     pDblY->IncreaseRef();
1326
1327     in.push_back(pDblT);
1328     in.push_back(pDblY);
1329
1330     for (int i = 0; i < (int)m_JacArgs.size(); i++)
1331     {
1332         m_JacArgs[i]->IncreaseRef();
1333         in.push_back(m_JacArgs[i]);
1334     }
1335
1336     try
1337     {
1338         // new std::string("") is delete in destructor of ast::CommentExp
1339         m_pCallJacFunction->invoke(in, opt, iRetCount, out, ast::CommentExp(Location(), new std::string("")));
1340     }
1341     catch (const ast::InternalError& ie)
1342     {
1343         for (int i = 0; i < (int)m_JacArgs.size(); i++)
1344         {
1345             m_JacArgs[i]->DecreaseRef();
1346         }
1347
1348         throw ie;
1349     }
1350
1351     for (int i = 0; i < (int)m_JacArgs.size(); i++)
1352     {
1353         m_JacArgs[i]->DecreaseRef();
1354     }
1355
1356     out[0]->IncreaseRef();
1357     pDblT->DecreaseRef();
1358     pDblY->DecreaseRef();
1359
1360     if (pDblT->isDeletable())
1361     {
1362         delete pDblT;
1363     }
1364
1365     if (pDblY->isDeletable())
1366     {
1367         delete pDblY;
1368     }
1369
1370     if (out.size() != iRetCount)
1371     {
1372         const char* pstrName = m_pCallJacFunction->getName().c_str();
1373         sprintf(errorMsg, _("%s: Wrong number of output argument(s): %d expected.\n"), pstrName, iRetCount);
1374         throw ast::InternalError(errorMsg);
1375     }
1376
1377     out[0]->DecreaseRef();
1378     if (out[0]->isDouble() == false)
1379     {
1380         const char* pstrName = m_pCallJacFunction->getName().c_str();
1381         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Real matrix expected.\n"), pstrName, 1);
1382         throw ast::InternalError(errorMsg);
1383     }
1384
1385
1386     types::Double* pDblOut = out[0]->getAs<types::Double>();
1387     int iSizeOut = pDblOut->getSize();
1388
1389     if (iSizeOut > iMaxSize)
1390     {
1391         const char* pstrName = m_pCallJacFunction->getName().c_str();
1392         sprintf(errorMsg, _("%s: Wrong size for output argument #%d: A size less or equal than %d expected.\n"), pstrName, 1, iMaxSize);
1393         throw ast::InternalError(errorMsg);
1394     }
1395
1396     C2F(dcopy)(&iSizeOut, pDblOut->get(), &one, J, &one);
1397 }
1398
1399 void DifferentialEquationFunctions::callMacroG(int* n, double* t, double* y, int* ng, double* gout)
1400 {
1401     char errorMsg[256];
1402     int iRetCount   = 1;
1403     int one         = 1;
1404
1405     types::typed_list in;
1406     types::typed_list out;
1407     types::optional_list opt;
1408
1409     types::Double* pDblY = new types::Double(m_odeYRows, m_odeYCols);
1410     pDblY->set(y);
1411     types::Double* pDblT = new types::Double(*t);
1412
1413     pDblT->IncreaseRef();
1414     pDblY->IncreaseRef();
1415
1416     in.push_back(pDblT);
1417     in.push_back(pDblY);
1418
1419     for (int i = 0; i < (int)m_odeGArgs.size(); i++)
1420     {
1421         m_odeGArgs[i]->IncreaseRef();
1422         in.push_back(m_odeGArgs[i]);
1423     }
1424
1425     try
1426     {
1427         // new std::string("") is delete in destructor of ast::CommentExp
1428         m_pCallGFunction->invoke(in, opt, iRetCount, out, ast::CommentExp(Location(), new std::string("")));
1429     }
1430     catch (const ast::InternalError& ie)
1431     {
1432         for (int i = 0; i < (int)m_odeGArgs.size(); i++)
1433         {
1434             m_odeGArgs[i]->DecreaseRef();
1435         }
1436
1437         throw ie;
1438     }
1439
1440     for (int i = 0; i < (int)m_odeGArgs.size(); i++)
1441     {
1442         m_odeGArgs[i]->DecreaseRef();
1443     }
1444
1445     if (out.size() != iRetCount)
1446     {
1447         const char* pstrName = m_pCallGFunction->getName().c_str();
1448         sprintf(errorMsg, _("%s: Wrong number of output argument(s): %d expected.\n"), pstrName, iRetCount);
1449         throw ast::InternalError(errorMsg);
1450     }
1451
1452     out[0]->IncreaseRef();
1453
1454     pDblT->DecreaseRef();
1455     pDblY->DecreaseRef();
1456
1457     if (pDblT->isDeletable())
1458     {
1459         delete pDblT;
1460     }
1461
1462     if (pDblY->isDeletable())
1463     {
1464         delete pDblY;
1465     }
1466
1467     out[0]->DecreaseRef();
1468     if (out[0]->isDouble() == false)
1469     {
1470         const char* pstrName = m_pCallGFunction->getName().c_str();
1471         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Real matrix expected.\n"), pstrName, 1);
1472         throw ast::InternalError(errorMsg);
1473     }
1474
1475     C2F(dcopy)(ng, out[0]->getAs<types::Double>()->get(), &one, gout, &one);
1476     if (out[0]->isDeletable())
1477     {
1478         delete out[0];
1479     }
1480 }
1481
1482 double DifferentialEquationFunctions::callIntgMacroF(double* t)
1483 {
1484     char errorMsg[256];
1485     int one         = 1;
1486     int iRetCount   = 1;
1487
1488     types::typed_list in;
1489     types::typed_list out;
1490     types::optional_list opt;
1491
1492     // create input args
1493     types::Double* pDblT = new types::Double(*t);
1494     pDblT->IncreaseRef();
1495
1496     // push_back
1497     in.push_back(pDblT);
1498
1499     for (int i = 0; i < (int)m_FArgs.size(); i++)
1500     {
1501         m_FArgs[i]->IncreaseRef();
1502         in.push_back(m_FArgs[i]);
1503     }
1504
1505     try
1506     {
1507         // new std::string("") is delete in destructor of ast::CommentExp
1508         m_pCallFFunction->invoke(in, opt, iRetCount, out, ast::CommentExp(Location(), new std::string("")));
1509     }
1510     catch (const ast::InternalError& ie)
1511     {
1512         for (int i = 0; i < (int)m_FArgs.size(); i++)
1513         {
1514             m_FArgs[i]->DecreaseRef();
1515         }
1516
1517         throw ie;
1518     }
1519
1520     for (int i = 0; i < (int)m_FArgs.size(); i++)
1521     {
1522         m_FArgs[i]->DecreaseRef();
1523     }
1524
1525     if (out.size() != iRetCount)
1526     {
1527         const char* pstrName = m_pCallFFunction->getName().c_str();
1528         sprintf(errorMsg, _("%s: Wrong number of output argument(s): %d expected.\n"), pstrName, iRetCount);
1529         throw ast::InternalError(errorMsg);
1530     }
1531
1532     out[0]->IncreaseRef();
1533
1534     pDblT->DecreaseRef();
1535     if (pDblT->isDeletable())
1536     {
1537         delete pDblT;
1538     }
1539
1540     out[0]->DecreaseRef();
1541     if (out[0]->isDouble() == false)
1542     {
1543         const char* pstrName = m_pCallFFunction->getName().c_str();
1544         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Real matrix expected.\n"), pstrName, 1);
1545         throw ast::InternalError(errorMsg);
1546
1547     }
1548
1549     types::Double* pDblOut = out[0]->getAs<types::Double>();
1550     if (pDblOut->getSize() != 1)
1551     {
1552         const char* pstrName = m_pCallFFunction->getName().c_str();
1553         sprintf(errorMsg, _("%s: Wrong size for output argument #%d: A Scalar expected.\n"), pstrName, 1);
1554         throw ast::InternalError(errorMsg);
1555     }
1556
1557     double res = pDblOut->get(0);
1558     if (out[0]->isDeletable())
1559     {
1560         delete out[0];
1561     }
1562
1563     return res;
1564 }
1565
1566 double DifferentialEquationFunctions::callInt2dMacroF(double* x, double* y)
1567 {
1568     char errorMsg[256];
1569     int one         = 1;
1570     int iRetCount   = 1;
1571
1572     types::typed_list in;
1573     types::typed_list out;
1574     types::optional_list opt;
1575
1576     // create input args
1577     types::Double* pDblX = new types::Double(*x);
1578     pDblX->IncreaseRef();
1579     types::Double* pDblY = new types::Double(*y);
1580     pDblY->IncreaseRef();
1581
1582     // push_back
1583     in.push_back(pDblX);
1584     in.push_back(pDblY);
1585
1586     for (int i = 0; i < (int)m_FArgs.size(); i++)
1587     {
1588         m_FArgs[i]->IncreaseRef();
1589         in.push_back(m_FArgs[i]);
1590     }
1591
1592     try
1593     {
1594         // new std::string("") is delete in destructor of ast::CommentExp
1595         m_pCallFFunction->invoke(in, opt, iRetCount, out, ast::CommentExp(Location(), new std::string("")));
1596     }
1597     catch (const ast::InternalError& ie)
1598     {
1599         for (int i = 0; i < (int)m_FArgs.size(); i++)
1600         {
1601             m_FArgs[i]->DecreaseRef();
1602         }
1603
1604         throw ie;
1605     }
1606
1607     for (int i = 0; i < (int)m_FArgs.size(); i++)
1608     {
1609         m_FArgs[i]->DecreaseRef();
1610     }
1611
1612     if (out.size() != iRetCount)
1613     {
1614         const char* pstrName = m_pCallFFunction->getName().c_str();
1615         sprintf(errorMsg, _("%s: Wrong number of output argument(s): %d expected.\n"), pstrName, iRetCount);
1616         throw ast::InternalError(errorMsg);
1617     }
1618
1619     out[0]->IncreaseRef();
1620
1621     pDblX->DecreaseRef();
1622     if (pDblX->isDeletable())
1623     {
1624         delete pDblX;
1625     }
1626
1627     pDblY->DecreaseRef();
1628     if (pDblY->isDeletable())
1629     {
1630         delete pDblY;
1631     }
1632
1633     out[0]->DecreaseRef();
1634     if (out[0]->isDouble() == false)
1635     {
1636         const char* pstrName = m_pCallFFunction->getName().c_str();
1637         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Real matrix expected.\n"), pstrName, 1);
1638         throw ast::InternalError(errorMsg);
1639     }
1640
1641     types::Double* pDblOut = out[0]->getAs<types::Double>();
1642     if (pDblOut->getSize() != 1)
1643     {
1644         const char* pstrName = m_pCallFFunction->getName().c_str();
1645         sprintf(errorMsg, _("%s: Wrong size for output argument #%d: A Scalar expected.\n"), pstrName, 1);
1646         throw ast::InternalError(errorMsg);
1647     }
1648
1649     double res = pDblOut->get(0);
1650     if (out[0]->isDeletable())
1651     {
1652         delete out[0];
1653     }
1654
1655     return res;
1656 }
1657
1658 void DifferentialEquationFunctions::callInt3dMacroF(double* xyz, int* numfun, double* funvls)
1659 {
1660     char errorMsg[256];
1661     int one         = 1;
1662     int iRetCount   = 1;
1663
1664     types::typed_list in;
1665     types::typed_list out;
1666     types::optional_list opt;
1667
1668     // create input args
1669     types::Double* pDblXYZ = new types::Double(3, 1);
1670     pDblXYZ->set(xyz);
1671     pDblXYZ->IncreaseRef();
1672     types::Double* pDblNumfun = new types::Double(*numfun);
1673     pDblNumfun->IncreaseRef();
1674
1675     // push_back
1676     in.push_back(pDblXYZ);
1677     in.push_back(pDblNumfun);
1678
1679     for (int i = 0; i < (int)m_FArgs.size(); i++)
1680     {
1681         m_FArgs[i]->IncreaseRef();
1682         in.push_back(m_FArgs[i]);
1683     }
1684
1685     try
1686     {
1687         // new std::string("") is delete in destructor of ast::CommentExp
1688         m_pCallFFunction->invoke(in, opt, iRetCount, out, ast::CommentExp(Location(), new std::string("")));
1689     }
1690     catch (const ast::InternalError& ie)
1691     {
1692         for (int i = 0; i < (int)m_FArgs.size(); i++)
1693         {
1694             m_FArgs[i]->DecreaseRef();
1695         }
1696
1697         throw ie;
1698     }
1699
1700     for (int i = 0; i < (int)m_FArgs.size(); i++)
1701     {
1702         m_FArgs[i]->DecreaseRef();
1703     }
1704
1705     if (out.size() != iRetCount)
1706     {
1707         const char* pstrName = m_pCallFFunction->getName().c_str();
1708         sprintf(errorMsg, _("%s: Wrong number of output argument(s): %d expected.\n"), pstrName, iRetCount);
1709         throw ast::InternalError(errorMsg);
1710     }
1711
1712     out[0]->IncreaseRef();
1713
1714     pDblXYZ->DecreaseRef();
1715     if (pDblXYZ->isDeletable())
1716     {
1717         delete pDblXYZ;
1718     }
1719
1720     pDblNumfun->DecreaseRef();
1721     if (pDblNumfun->isDeletable())
1722     {
1723         delete pDblNumfun;
1724     }
1725
1726     out[0]->DecreaseRef();
1727     if (out[0]->isDouble() == false)
1728     {
1729         const char* pstrName = m_pCallFFunction->getName().c_str();
1730         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Real matrix expected.\n"), pstrName, 1);
1731         throw ast::InternalError(errorMsg);
1732     }
1733
1734     types::Double* pDblOut = out[0]->getAs<types::Double>();
1735     if (pDblOut->getSize() != *numfun)
1736     {
1737         const char* pstrName = m_pCallFFunction->getName().c_str();
1738         sprintf(errorMsg, _("%s: Wrong size for output argument #%d: Matrix of size %d expected.\n"), pstrName, 1, *numfun);
1739         throw ast::InternalError(errorMsg);
1740     }
1741
1742     C2F(dcopy)(numfun, pDblOut->get(), &one, funvls, &one);
1743     if (out[0]->isDeletable())
1744     {
1745         delete out[0];
1746     }
1747 }
1748
1749 void DifferentialEquationFunctions::callFevalMacroF(int* nn, double* x1, double* x2, double* xres, int* itype)
1750 {
1751     char errorMsg[256];
1752     int one         = 1;
1753     int iRetCount   = 1;
1754
1755     types::typed_list in;
1756     types::typed_list out;
1757     types::optional_list opt;
1758
1759     types::Double* pDblX = NULL;
1760     types::Double* pDblY = NULL;
1761
1762     // create input args
1763
1764     pDblX = new types::Double(*x1);
1765     pDblX->IncreaseRef();
1766     in.push_back(pDblX);
1767
1768     if (*nn == 2)
1769     {
1770         pDblY = new types::Double(*x2);
1771         pDblY->IncreaseRef();
1772         in.push_back(pDblY);
1773     }
1774
1775     for (int i = 0; i < (int)m_FArgs.size(); i++)
1776     {
1777         m_FArgs[i]->IncreaseRef();
1778         in.push_back(m_FArgs[i]);
1779     }
1780     try
1781     {
1782         // new std::string("") is delete in destructor of ast::CommentExp
1783         m_pCallFFunction->invoke(in, opt, iRetCount, out, ast::CommentExp(Location(), new std::string("")));
1784     }
1785     catch (const ast::InternalError& ie)
1786     {
1787         for (int i = 0; i < (int)m_FArgs.size(); i++)
1788         {
1789             m_FArgs[i]->DecreaseRef();
1790         }
1791
1792         throw ie;
1793     }
1794
1795     for (int i = 0; i < (int)m_FArgs.size(); i++)
1796     {
1797         m_FArgs[i]->DecreaseRef();
1798     }
1799
1800     if (out.size() != iRetCount)
1801     {
1802         const char* pstrName = m_pCallFFunction->getName().c_str();
1803         sprintf(errorMsg, _("%s: Wrong number of output argument(s): %d expected.\n"), pstrName, iRetCount);
1804         throw ast::InternalError(errorMsg);
1805     }
1806
1807     out[0]->IncreaseRef();
1808
1809     pDblX->DecreaseRef();
1810     if (pDblX->isDeletable())
1811     {
1812         delete pDblX;
1813     }
1814
1815     if (*nn == 2)
1816     {
1817         pDblY->DecreaseRef();
1818         if (pDblY->isDeletable())
1819         {
1820             delete pDblY;
1821         }
1822     }
1823
1824     out[0]->DecreaseRef();
1825     if (out[0]->isDouble() == false)
1826     {
1827         const char* pstrName = m_pCallFFunction->getName().c_str();
1828         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Real matrix expected.\n"), pstrName, 1);
1829         throw ast::InternalError(errorMsg);
1830
1831     }
1832
1833     types::Double* pDblOut = out[0]->getAs<types::Double>();
1834     if (pDblOut->getSize() != 1)
1835     {
1836         const char* pstrName = m_pCallFFunction->getName().c_str();
1837         sprintf(errorMsg, _("%s: Wrong size for output argument #%d: A Scalar expected.\n"), pstrName, 1);
1838         throw ast::InternalError(errorMsg);
1839     }
1840
1841     if (pDblOut->isComplex())
1842     {
1843         *itype = 1;
1844         xres[0] = pDblOut->get(0);
1845         xres[1] = pDblOut->getImg(0);
1846     }
1847     else
1848     {
1849         *itype = 0;
1850         xres[0] = pDblOut->get(0);
1851     }
1852
1853     if (out[0]->isDeletable())
1854     {
1855         delete out[0];
1856     }
1857 }
1858
1859 void DifferentialEquationFunctions::callBvodeMacroGsub(int* i, double* z, double* g)
1860 {
1861     char errorMsg[256];
1862     int one         = 1;
1863     int iRetCount   = 1;
1864
1865     types::typed_list in;
1866     types::typed_list out;
1867     types::optional_list opt;
1868
1869     types::Double* pDblI = NULL;
1870     types::Double* pDblZ = NULL;
1871
1872     pDblI = new types::Double(*i);
1873     pDblI->IncreaseRef();
1874     in.push_back(pDblI);
1875
1876     pDblZ = new types::Double(m_bvodeM, 1);
1877     pDblZ->set(z);
1878     pDblZ->IncreaseRef();
1879     in.push_back(pDblZ);
1880
1881
1882     for (int i = 0; i < (int)m_GsubArgs.size(); i++)
1883     {
1884         m_GsubArgs[i]->IncreaseRef();
1885         in.push_back(m_GsubArgs[i]);
1886     }
1887
1888     try
1889     {
1890         // new std::string("") is delete in destructor of ast::CommentExp
1891         m_pCallGsubFunction->invoke(in, opt, iRetCount, out, ast::CommentExp(Location(), new std::string("")));
1892     }
1893     catch (const ast::InternalError& ie)
1894     {
1895         for (int i = 0; i < (int)m_GsubArgs.size(); i++)
1896         {
1897             m_GsubArgs[i]->DecreaseRef();
1898         }
1899
1900         throw ie;
1901     }
1902
1903     for (int i = 0; i < (int)m_GsubArgs.size(); i++)
1904     {
1905         m_GsubArgs[i]->DecreaseRef();
1906     }
1907
1908     if (out.size() != iRetCount)
1909     {
1910         const char* pstrName = m_pCallGsubFunction->getName().c_str();
1911         sprintf(errorMsg, _("%s: Wrong number of output argument(s): %d expected.\n"), pstrName, iRetCount);
1912         throw ast::InternalError(errorMsg);
1913     }
1914
1915     out[0]->IncreaseRef();
1916
1917     pDblI->DecreaseRef();
1918     if (pDblI->isDeletable())
1919     {
1920         delete pDblI;
1921     }
1922
1923     pDblZ->DecreaseRef();
1924     if (pDblZ->isDeletable())
1925     {
1926         delete pDblZ;
1927     }
1928
1929     out[0]->DecreaseRef();
1930     if (out[0]->isDouble() == false)
1931     {
1932         const char* pstrName = m_pCallGsubFunction->getName().c_str();
1933         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Real matrix expected.\n"), pstrName, 1);
1934         throw ast::InternalError(errorMsg);
1935     }
1936
1937     types::Double* pDblOut = out[0]->getAs<types::Double>();
1938     if (pDblOut->getSize() != 1)
1939     {
1940         const char* pstrName = m_pCallGsubFunction->getName().c_str();
1941         sprintf(errorMsg, _("%s: Wrong size for output argument #%d: A Scalar expected.\n"), pstrName, 1);
1942         throw ast::InternalError(errorMsg);
1943     }
1944
1945     *g = pDblOut->get(0);
1946     if (out[0]->isDeletable())
1947     {
1948         delete out[0];
1949     }
1950 }
1951
1952 void DifferentialEquationFunctions::callBvodeMacroDgsub(int* i, double* z, double* g)
1953 {
1954     char errorMsg[256];
1955     int one         = 1;
1956     int iRetCount   = 1;
1957
1958     types::typed_list in;
1959     types::typed_list out;
1960     types::optional_list opt;
1961
1962     types::Double* pDblI = NULL;
1963     types::Double* pDblZ = NULL;
1964
1965     pDblI = new types::Double(*i);
1966     pDblI->IncreaseRef();
1967     in.push_back(pDblI);
1968
1969     pDblZ = new types::Double(m_bvodeM, 1);
1970     pDblZ->set(z);
1971     pDblZ->IncreaseRef();
1972     in.push_back(pDblZ);
1973
1974     for (int i = 0; i < (int)m_DgsubArgs.size(); i++)
1975     {
1976         m_DgsubArgs[i]->IncreaseRef();
1977         in.push_back(m_DgsubArgs[i]);
1978     }
1979
1980     try
1981     {
1982         // new std::string("") is delete in destructor of ast::CommentExp
1983         m_pCallDgsubFunction->invoke(in, opt, iRetCount, out, ast::CommentExp(Location(), new std::string("")));
1984     }
1985     catch (const ast::InternalError& ie)
1986     {
1987         for (int i = 0; i < (int)m_DgsubArgs.size(); i++)
1988         {
1989             m_DgsubArgs[i]->DecreaseRef();
1990         }
1991
1992         throw ie;
1993     }
1994
1995     for (int i = 0; i < (int)m_DgsubArgs.size(); i++)
1996     {
1997         m_DgsubArgs[i]->DecreaseRef();
1998     }
1999
2000     if (out.size() != iRetCount)
2001     {
2002         const char* pstrName = m_pCallDgsubFunction->getName().c_str();
2003         sprintf(errorMsg, _("%s: Wrong number of output argument(s): %d expected.\n"), pstrName, iRetCount);
2004         throw ast::InternalError(errorMsg);
2005     }
2006
2007     out[0]->IncreaseRef();
2008
2009     pDblI->DecreaseRef();
2010     if (pDblI->isDeletable())
2011     {
2012         delete pDblI;
2013     }
2014
2015     pDblZ->DecreaseRef();
2016     if (pDblZ->isDeletable())
2017     {
2018         delete pDblZ;
2019     }
2020
2021     out[0]->DecreaseRef();
2022     if (out[0]->isDouble() == false)
2023     {
2024         const char* pstrName = m_pCallDgsubFunction->getName().c_str();
2025         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Real matrix expected.\n"), pstrName, 1);
2026         throw ast::InternalError(errorMsg);
2027     }
2028
2029     types::Double* pDblOut = out[0]->getAs<types::Double>();
2030     if (pDblOut->getSize() != m_bvodeM)
2031     {
2032         const char* pstrName = m_pCallDgsubFunction->getName().c_str();
2033         sprintf(errorMsg, _("%s: Wrong size for output argument #%d: A Matrix of size %d expected.\n"), pstrName, 1, m_bvodeM);
2034         throw ast::InternalError(errorMsg);
2035     }
2036
2037     C2F(dcopy)(&m_bvodeM, pDblOut->get(), &one, g, &one);
2038     if (out[0]->isDeletable())
2039     {
2040         delete out[0];
2041     }
2042 }
2043
2044 void DifferentialEquationFunctions::callBvodeMacroFsub(double* x, double* z, double* d)
2045 {
2046     char errorMsg[256];
2047     int one         = 1;
2048     int iRetCount   = 1;
2049
2050     types::typed_list in;
2051     types::typed_list out;
2052     types::optional_list opt;
2053
2054     types::Double* pDblX = NULL;
2055     types::Double* pDblZ = NULL;
2056
2057     pDblX = new types::Double(*x);
2058     pDblX->IncreaseRef();
2059     in.push_back(pDblX);
2060
2061     pDblZ = new types::Double(m_bvodeM, 1);
2062     pDblZ->set(z);
2063     pDblZ->IncreaseRef();
2064     in.push_back(pDblZ);
2065
2066     for (int i = 0; i < (int)m_FsubArgs.size(); i++)
2067     {
2068         m_FsubArgs[i]->IncreaseRef();
2069         in.push_back(m_FsubArgs[i]);
2070     }
2071
2072     try
2073     {
2074         // new std::string("") is delete in destructor of ast::CommentExp
2075         m_pCallFsubFunction->invoke(in, opt, iRetCount, out, ast::CommentExp(Location(), new std::string("")));
2076     }
2077     catch (const ast::InternalError& ie)
2078     {
2079         for (int i = 0; i < (int)m_FsubArgs.size(); i++)
2080         {
2081             m_FsubArgs[i]->DecreaseRef();
2082         }
2083
2084         throw ie;
2085     }
2086
2087     for (int i = 0; i < (int)m_FsubArgs.size(); i++)
2088     {
2089         m_FsubArgs[i]->DecreaseRef();
2090     }
2091
2092     if (out.size() != iRetCount)
2093     {
2094         const char* pstrName = m_pCallFsubFunction->getName().c_str();
2095         sprintf(errorMsg, _("%s: Wrong number of output argument(s): %d expected.\n"), pstrName, iRetCount);
2096         throw ast::InternalError(errorMsg);
2097     }
2098
2099     out[0]->IncreaseRef();
2100
2101     pDblX->DecreaseRef();
2102     if (pDblX->isDeletable())
2103     {
2104         delete pDblX;
2105     }
2106
2107     pDblZ->DecreaseRef();
2108     if (pDblZ->isDeletable())
2109     {
2110         delete pDblZ;
2111     }
2112
2113     out[0]->DecreaseRef();
2114     if (out[0]->isDouble() == false)
2115     {
2116         const char* pstrName = m_pCallFsubFunction->getName().c_str();
2117         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Real matrix expected.\n"), pstrName, 1);
2118         throw ast::InternalError(errorMsg);
2119     }
2120
2121     types::Double* pDblOut = out[0]->getAs<types::Double>();
2122     if (pDblOut->getSize() != m_bvodeN)
2123     {
2124         const char* pstrName = m_pCallFsubFunction->getName().c_str();
2125         sprintf(errorMsg, _("%s: Wrong size for output argument #%d: A matrix of size %d expected.\n"), pstrName, 1, m_bvodeN);
2126         throw ast::InternalError(errorMsg);
2127     }
2128
2129     C2F(dcopy)(&m_bvodeN, pDblOut->get(), &one, d, &one);
2130     if (out[0]->isDeletable())
2131     {
2132         delete out[0];
2133     }
2134 }
2135
2136 void DifferentialEquationFunctions::callBvodeMacroDfsub(double* x, double* z, double* d)
2137 {
2138     char errorMsg[256];
2139     int one         = 1;
2140     int iRetCount   = 1;
2141
2142     types::typed_list in;
2143     types::typed_list out;
2144     types::optional_list opt;
2145
2146     types::Double* pDblX = NULL;
2147     types::Double* pDblZ = NULL;
2148
2149     pDblX = new types::Double(*x);
2150     pDblX->IncreaseRef();
2151     in.push_back(pDblX);
2152
2153     pDblZ = new types::Double(m_bvodeM, 1);
2154     pDblZ->set(z);
2155     pDblZ->IncreaseRef();
2156     in.push_back(pDblZ);
2157
2158     for (int i = 0; i < (int)m_DfsubArgs.size(); i++)
2159     {
2160         m_DfsubArgs[i]->IncreaseRef();
2161         in.push_back(m_DfsubArgs[i]);
2162     }
2163
2164     try
2165     {
2166         // new std::string("") is delete in destructor of ast::CommentExp
2167         m_pCallDfsubFunction->invoke(in, opt, iRetCount, out, ast::CommentExp(Location(), new std::string("")));
2168     }
2169     catch (const ast::InternalError& ie)
2170     {
2171         for (int i = 0; i < (int)m_DfsubArgs.size(); i++)
2172         {
2173             m_DfsubArgs[i]->DecreaseRef();
2174         }
2175
2176         throw ie;
2177     }
2178
2179     for (int i = 0; i < (int)m_DfsubArgs.size(); i++)
2180     {
2181         m_DfsubArgs[i]->DecreaseRef();
2182     }
2183
2184     if (out.size() != iRetCount)
2185     {
2186         const char* pstrName = m_pCallDfsubFunction->getName().c_str();
2187         sprintf(errorMsg, _("%s: Wrong number of output argument(s): %d expected.\n"), pstrName, iRetCount);
2188         throw ast::InternalError(errorMsg);
2189     }
2190
2191     out[0]->IncreaseRef();
2192
2193     pDblX->DecreaseRef();
2194     if (pDblX->isDeletable())
2195     {
2196         delete pDblX;
2197     }
2198
2199     pDblZ->DecreaseRef();
2200     if (pDblZ->isDeletable())
2201     {
2202         delete pDblZ;
2203     }
2204
2205     out[0]->DecreaseRef();
2206     if (out[0]->isDouble() == false)
2207     {
2208         const char* pstrName = m_pCallDfsubFunction->getName().c_str();
2209         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Real matrix expected.\n"), pstrName, 1);
2210         throw ast::InternalError(errorMsg);
2211     }
2212
2213     types::Double* pDblOut = out[0]->getAs<types::Double>();
2214     int size = m_bvodeN * m_bvodeM;
2215     if (pDblOut->getSize() != size)
2216     {
2217         const char* pstrName = m_pCallDfsubFunction->getName().c_str();
2218         sprintf(errorMsg, _("%s: Wrong size for output argument #%d: A matrix of size %d expected.\n"), pstrName, 1, size);
2219         throw ast::InternalError(errorMsg);
2220     }
2221
2222     C2F(dcopy)(&size, pDblOut->get(), &one, d, &one);
2223     if (out[0]->isDeletable())
2224     {
2225         delete out[0];
2226     }
2227 }
2228
2229 void DifferentialEquationFunctions::callBvodeMacroGuess(double* x, double* z, double* d)
2230 {
2231     char errorMsg[256];
2232     int one         = 1;
2233     int iRetCount   = 2;
2234
2235     types::typed_list in;
2236     types::typed_list out;
2237     types::optional_list opt;
2238
2239     types::Double* pDblX = NULL;
2240
2241     pDblX = new types::Double(*x);
2242     pDblX->IncreaseRef();
2243     in.push_back(pDblX);
2244
2245     for (int i = 0; i < (int)m_GuessArgs.size(); i++)
2246     {
2247         m_GuessArgs[i]->IncreaseRef();
2248         in.push_back(m_GuessArgs[i]);
2249     }
2250
2251     try
2252     {
2253         // new std::string("") is delete in destructor of ast::CommentExp
2254         m_pCallGuessFunction->invoke(in, opt, iRetCount, out, ast::CommentExp(Location(), new std::string("")));
2255     }
2256     catch (const ast::InternalError& ie)
2257     {
2258         for (int i = 0; i < (int)m_GuessArgs.size(); i++)
2259         {
2260             m_GuessArgs[i]->DecreaseRef();
2261         }
2262
2263         throw ie;
2264     }
2265
2266     for (int i = 0; i < (int)m_GuessArgs.size(); i++)
2267     {
2268         m_GuessArgs[i]->DecreaseRef();
2269     }
2270
2271     if (out.size() != iRetCount)
2272     {
2273         const char* pstrName = m_pCallGuessFunction->getName().c_str();
2274         sprintf(errorMsg, _("%s: Wrong number of output argument(s): %d expected.\n"), pstrName, iRetCount);
2275         throw ast::InternalError(errorMsg);
2276     }
2277
2278     out[0]->IncreaseRef();
2279     out[1]->IncreaseRef();
2280
2281     pDblX->DecreaseRef();
2282     if (pDblX->isDeletable())
2283     {
2284         delete pDblX;
2285     }
2286
2287     out[0]->DecreaseRef();
2288     if (out[0]->isDouble() == false)
2289     {
2290         const char* pstrName = m_pCallGuessFunction->getName().c_str();
2291         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Real matrix expected.\n"), pstrName, 1);
2292         throw ast::InternalError(errorMsg);
2293     }
2294
2295     out[1]->DecreaseRef();
2296     if (out[1]->isDouble() == false)
2297     {
2298         const char* pstrName = m_pCallGuessFunction->getName().c_str();
2299         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Real matrix expected.\n"), pstrName, 2);
2300         throw ast::InternalError(errorMsg);
2301     }
2302
2303     types::Double* pDblOutZ = out[0]->getAs<types::Double>();
2304     if (pDblOutZ->getSize() != m_bvodeM)
2305     {
2306         const char* pstrName = m_pCallGuessFunction->getName().c_str();
2307         sprintf(errorMsg, _("%s: Wrong size for output argument #%d: A matrix of size %d expected.\n"), pstrName, 1, m_bvodeM);
2308         throw ast::InternalError(errorMsg);
2309     }
2310
2311     types::Double* pDblOutD = out[1]->getAs<types::Double>();
2312     if (pDblOutD->getSize() != m_bvodeN)
2313     {
2314         const char* pstrName = m_pCallGuessFunction->getName().c_str();
2315         sprintf(errorMsg, _("%s: Wrong size for output argument #%d: A matrix of size %d expected.\n"), pstrName, 1, m_bvodeN);
2316         throw ast::InternalError(errorMsg);
2317     }
2318
2319     C2F(dcopy)(&m_bvodeM, pDblOutZ->get(), &one, z, &one);
2320     C2F(dcopy)(&m_bvodeN, pDblOutD->get(), &one, d, &one);
2321     if (out[0]->isDeletable())
2322     {
2323         delete out[0];
2324     }
2325
2326     if (out[1]->isDeletable())
2327     {
2328         delete out[1];
2329     }
2330 }
2331
2332 void DifferentialEquationFunctions::callImplMacroF(int* neq, double* t, double* y, double*s, double* r, int* ires)
2333 {
2334     char errorMsg[256];
2335     int one         = 1;
2336     int iRetCount   = 1;
2337
2338     *ires = 2;
2339
2340     types::typed_list in;
2341     types::typed_list out;
2342     types::optional_list opt;
2343
2344     types::Double* pDblT = new types::Double(*t);
2345     pDblT->IncreaseRef();
2346     in.push_back(pDblT);
2347
2348     types::Double* pDblY = new types::Double(*neq, 1);
2349     pDblY->set(y);
2350     pDblY->IncreaseRef();
2351     in.push_back(pDblY);
2352
2353     types::Double* pDblS = new types::Double(*neq, 1);
2354     pDblS->set(s);
2355     pDblS->IncreaseRef();
2356     in.push_back(pDblS);
2357
2358     for (int i = 0; i < (int)m_FArgs.size(); i++)
2359     {
2360         m_FArgs[i]->IncreaseRef();
2361         in.push_back(m_FArgs[i]);
2362     }
2363
2364     try
2365     {
2366         // new std::string("") is delete in destructor of ast::CommentExp
2367         m_pCallFFunction->invoke(in, opt, iRetCount, out, ast::CommentExp(Location(), new std::string("")));
2368     }
2369     catch (const ast::InternalError& ie)
2370     {
2371         for (int i = 0; i < (int)m_FArgs.size(); i++)
2372         {
2373             m_FArgs[i]->DecreaseRef();
2374         }
2375
2376         throw ie;
2377     }
2378
2379     if (out.size() != iRetCount)
2380     {
2381         const char* pstrName = m_pCallFFunction->getName().c_str();
2382         sprintf(errorMsg, _("%s: Wrong number of output argument(s): %d expected.\n"), pstrName, iRetCount);
2383         throw ast::InternalError(errorMsg);
2384     }
2385
2386     out[0]->IncreaseRef();
2387
2388     pDblT->DecreaseRef();
2389     if (pDblT->isDeletable())
2390     {
2391         delete pDblT;
2392     }
2393
2394     pDblY->DecreaseRef();
2395     if (pDblY->isDeletable())
2396     {
2397         delete pDblY;
2398     }
2399
2400     pDblS->DecreaseRef();
2401     if (pDblS->isDeletable())
2402     {
2403         delete pDblS;
2404     }
2405
2406     out[0]->DecreaseRef();
2407     if (out[0]->isDouble() == false)
2408     {
2409         const char* pstrName = m_pCallFFunction->getName().c_str();
2410         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Matrix expected.\n"), pstrName, 1);
2411         throw ast::InternalError(errorMsg);
2412     }
2413
2414     types::Double* pDblOutR = out[0]->getAs<types::Double>();
2415     if (pDblOutR->getSize() != *neq)
2416     {
2417         const char* pstrName = m_pCallFFunction->getName().c_str();
2418         sprintf(errorMsg, _("%s: Wrong size for output argument #%d: A Matrix of size %d expected.\n"), pstrName, 1, *neq);
2419         throw ast::InternalError(errorMsg);
2420     }
2421
2422     C2F(dcopy)(neq, pDblOutR->get(), &one, r, &one);
2423     *ires = 1;
2424     if (out[0]->isDeletable())
2425     {
2426         delete out[0];
2427     }
2428 }
2429
2430 void DifferentialEquationFunctions::callImplMacroG(int* neq, double* t, double* y, double* ml, double* mu, double* p, int* nrowp)
2431 {
2432     char errorMsg[256];
2433     int one         = 1;
2434     int iRetCount   = 1;
2435
2436     types::typed_list in;
2437     types::typed_list out;
2438     types::optional_list opt;
2439
2440     types::Double* pDblT = new types::Double(*t);
2441     pDblT->IncreaseRef();
2442     in.push_back(pDblT);
2443
2444     types::Double* pDblY = new types::Double(*neq, 1);
2445     pDblY->set(y);
2446     pDblY->IncreaseRef();
2447     in.push_back(pDblY);
2448
2449     types::Double* pDblP = new types::Double(*nrowp, *neq);
2450     pDblP->set(p);
2451     pDblP->IncreaseRef();
2452     in.push_back(pDblP);
2453
2454     for (int i = 0; i < (int)m_odeGArgs.size(); i++)
2455     {
2456         m_odeGArgs[i]->IncreaseRef();
2457         in.push_back(m_odeGArgs[i]);
2458     }
2459
2460     try
2461     {
2462         // new std::string("") is delete in destructor of ast::CommentExp
2463         m_pCallGFunction->invoke(in, opt, iRetCount, out, ast::CommentExp(Location(), new std::string("")));
2464     }
2465     catch (const ast::InternalError& ie)
2466     {
2467         for (int i = 0; i < (int)m_odeGArgs.size(); i++)
2468         {
2469             m_odeGArgs[i]->DecreaseRef();
2470         }
2471
2472         throw ie;
2473     }
2474
2475     for (int i = 0; i < (int)m_odeGArgs.size(); i++)
2476     {
2477         m_odeGArgs[i]->DecreaseRef();
2478     }
2479
2480     if (out.size() != iRetCount)
2481     {
2482         const char* pstrName = m_pCallGFunction->getName().c_str();
2483         sprintf(errorMsg, _("%s: Wrong number of output argument(s): %d expected.\n"), pstrName, iRetCount);
2484         throw ast::InternalError(errorMsg);
2485     }
2486
2487     out[0]->IncreaseRef();
2488
2489     pDblT->DecreaseRef();
2490     if (pDblT->isDeletable())
2491     {
2492         delete pDblT;
2493     }
2494
2495     pDblY->DecreaseRef();
2496     if (pDblY->isDeletable())
2497     {
2498         delete pDblY;
2499     }
2500
2501     out[0]->DecreaseRef();
2502     if (out[0]->isDouble() == false)
2503     {
2504         const char* pstrName = m_pCallGFunction->getName().c_str();
2505         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Real matrix expected.\n"), pstrName, 1);
2506         throw ast::InternalError(errorMsg);
2507     }
2508
2509     types::Double* pDblOutP = out[0]->getAs<types::Double>();
2510     if (pDblOutP->getCols() != *neq || pDblOutP->getRows() != *nrowp)
2511     {
2512         const char* pstrName = m_pCallGFunction->getName().c_str();
2513         sprintf(errorMsg, _("%s: Wrong size for output argument #%d: A matrix of size %d x %d expected.\n"), pstrName, 1, *neq, *nrowp);
2514         throw ast::InternalError(errorMsg);
2515     }
2516
2517     int size = *neq **nrowp;
2518     C2F(dcopy)(&size, pDblOutP->get(), &one, p, &one);
2519     if (out[0]->isDeletable())
2520     {
2521         delete out[0];
2522     }
2523 }
2524
2525 void DifferentialEquationFunctions::callImplMacroJac(int* neq, double* t, double* y, double* s, double* ml, double* mu, double* p, int* nrowp)
2526 {
2527     char errorMsg[256];
2528     int one         = 1;
2529     int iRetCount   = 1;
2530
2531     types::typed_list in;
2532     types::typed_list out;
2533     types::optional_list opt;
2534
2535     types::Double* pDblT = new types::Double(*t);
2536     pDblT->IncreaseRef();
2537     in.push_back(pDblT);
2538
2539     types::Double* pDblY = new types::Double(*neq, 1);
2540     pDblY->set(y);
2541     pDblY->IncreaseRef();
2542     in.push_back(pDblY);
2543
2544     types::Double* pDblS = new types::Double(*neq, 1);
2545     pDblS->set(s);
2546     pDblS->IncreaseRef();
2547     in.push_back(pDblS);
2548
2549     for (int i = 0; i < (int)m_JacArgs.size(); i++)
2550     {
2551         m_JacArgs[i]->IncreaseRef();
2552         in.push_back(m_JacArgs[i]);
2553     }
2554
2555     try
2556     {
2557         // new std::string("") is delete in destructor of ast::CommentExp
2558         m_pCallJacFunction->invoke(in, opt, iRetCount, out, ast::CommentExp(Location(), new std::string("")));
2559     }
2560     catch (const ast::InternalError& ie)
2561     {
2562         for (int i = 0; i < (int)m_JacArgs.size(); i++)
2563         {
2564             m_JacArgs[i]->DecreaseRef();
2565         }
2566
2567         throw ie;
2568     }
2569
2570     for (int i = 0; i < (int)m_JacArgs.size(); i++)
2571     {
2572         m_JacArgs[i]->DecreaseRef();
2573     }
2574
2575     if (out.size() != iRetCount)
2576     {
2577         const char* pstrName = m_pCallJacFunction->getName().c_str();
2578         sprintf(errorMsg, _("%s: Wrong number of output argument(s): %d expected.\n"), pstrName, iRetCount);
2579         throw ast::InternalError(errorMsg);
2580     }
2581
2582     out[0]->IncreaseRef();
2583
2584     pDblT->DecreaseRef();
2585     if (pDblT->isDeletable())
2586     {
2587         delete pDblT;
2588     }
2589
2590     pDblY->DecreaseRef();
2591     if (pDblY->isDeletable())
2592     {
2593         delete pDblY;
2594     }
2595
2596     pDblS->DecreaseRef();
2597     if (pDblS->isDeletable())
2598     {
2599         delete pDblS;
2600     }
2601
2602     out[0]->DecreaseRef();
2603     if (out[0]->isDouble() == false)
2604     {
2605         const char* pstrName = m_pCallJacFunction->getName().c_str();
2606         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Real matrix expected.\n"), pstrName, 1);
2607         throw ast::InternalError(errorMsg);
2608     }
2609
2610     types::Double* pDblOutP = out[0]->getAs<types::Double>();
2611     if (pDblOutP->getCols() != *neq || pDblOutP->getRows() != *nrowp)
2612     {
2613         const char* pstrName = m_pCallJacFunction->getName().c_str();
2614         sprintf(errorMsg, _("%s: Wrong size for output argument #%d: A matrix of size %d x %d expected.\n"), pstrName, 1, *neq, *nrowp);
2615         throw ast::InternalError(errorMsg);
2616     }
2617
2618     int size = *neq **nrowp;
2619     C2F(dcopy)(&size, pDblOutP->get(), &one, p, &one);
2620     if (out[0]->isDeletable())
2621     {
2622         delete out[0];
2623     }
2624 }
2625
2626 void DifferentialEquationFunctions::callDasslMacroF(double* t, double* y, double* ydot, double* delta, int* ires, double* rpar, int* ipar)
2627 {
2628     char errorMsg[256];
2629     int one         = 1;
2630     int iRetCount   = 2;
2631
2632     types::typed_list in;
2633     types::typed_list out;
2634     types::optional_list opt;
2635
2636     types::Double* pDblT = new types::Double(*t);
2637     pDblT->IncreaseRef();
2638     in.push_back(pDblT);
2639
2640     types::Double* pDblY = new types::Double(m_odeYRows, 1);
2641     pDblY->set(y);
2642     pDblY->IncreaseRef();
2643     in.push_back(pDblY);
2644
2645     types::Double* pDblYdot = new types::Double(m_odeYRows, 1);
2646     pDblYdot->set(ydot);
2647     pDblYdot->IncreaseRef();
2648     in.push_back(pDblYdot);
2649
2650     for (int i = 0; i < (int)m_FArgs.size(); i++)
2651     {
2652         m_FArgs[i]->IncreaseRef();
2653         in.push_back(m_FArgs[i]);
2654     }
2655
2656     try
2657     {
2658         // new std::string("") is delete in destructor of ast::CommentExp
2659         m_pCallFFunction->invoke(in, opt, iRetCount, out, ast::CommentExp(Location(), new std::string("")));
2660     }
2661     catch (const ast::InternalError& ie)
2662     {
2663         for (int i = 0; i < (int)m_FArgs.size(); i++)
2664         {
2665             m_FArgs[i]->DecreaseRef();
2666         }
2667
2668         throw ie;
2669     }
2670
2671     for (int i = 0; i < (int)m_FArgs.size(); i++)
2672     {
2673         m_FArgs[i]->DecreaseRef();
2674     }
2675
2676     if (out.size() != iRetCount)
2677     {
2678         const char* pstrName = m_pCallFFunction->getName().c_str();
2679         sprintf(errorMsg, _("%s: Wrong number of output argument(s): %d expected.\n"), pstrName, iRetCount);
2680         throw ast::InternalError(errorMsg);
2681     }
2682
2683     out[0]->IncreaseRef();
2684     out[1]->IncreaseRef();
2685
2686     pDblT->DecreaseRef();
2687     if (pDblT->isDeletable())
2688     {
2689         delete pDblT;
2690     }
2691
2692     pDblY->DecreaseRef();
2693     if (pDblY->isDeletable())
2694     {
2695         delete pDblY;
2696     }
2697
2698     pDblYdot->DecreaseRef();
2699     if (pDblYdot->isDeletable())
2700     {
2701         delete pDblYdot;
2702     }
2703
2704     out[0]->DecreaseRef();
2705     if (out[0]->isDouble() == false)
2706     {
2707         const char* pstrName = m_pCallFFunction->getName().c_str();
2708         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Real matrix expected.\n"), pstrName, 1);
2709         throw ast::InternalError(errorMsg);
2710     }
2711
2712     out[1]->DecreaseRef();
2713     if (out[1]->isDouble() == false)
2714     {
2715         const char* pstrName = m_pCallFFunction->getName().c_str();
2716         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Real matrix expected.\n"), pstrName, 2);
2717         throw ast::InternalError(errorMsg);
2718     }
2719
2720     types::Double* pDblOutDelta = out[0]->getAs<types::Double>();
2721     if (pDblOutDelta->getSize() != m_odeYRows)
2722     {
2723         const char* pstrName = m_pCallFFunction->getName().c_str();
2724         sprintf(errorMsg, _("%s: Wrong size for output argument #%d: A matrix of size %d expected.\n"), pstrName, 1, m_odeYRows);
2725         throw ast::InternalError(errorMsg);
2726     }
2727
2728     types::Double* pDblOutIres = out[1]->getAs<types::Double>();
2729     if (pDblOutIres->getSize() != 1)
2730     {
2731         const char* pstrName = m_pCallFFunction->getName().c_str();
2732         sprintf(errorMsg, _("%s: Wrong size for output argument #%d: A Scalar expected.\n"), pstrName, 2);
2733         throw ast::InternalError(errorMsg);
2734     }
2735
2736     C2F(dcopy)(&m_odeYRows, pDblOutDelta->get(), &one, delta, &one);
2737     *ires = (int)pDblOutIres->get(0);
2738
2739     if (out[0]->isDeletable())
2740     {
2741         delete out[0];
2742     }
2743
2744     if (out[1]->isDeletable())
2745     {
2746         delete out[1];
2747     }
2748 }
2749
2750 void DifferentialEquationFunctions::callDasslMacroJac(double* t, double* y, double* ydot, double* pd, double* cj, double* rpar, int* ipar)
2751 {
2752     char errorMsg[256];
2753     int one         = 1;
2754     int iRetCount   = 1;
2755
2756     types::typed_list in;
2757     types::typed_list out;
2758     types::optional_list opt;
2759
2760     types::Double* pDblT = new types::Double(*t);
2761     pDblT->IncreaseRef();
2762     in.push_back(pDblT);
2763
2764     types::Double* pDblY = new types::Double(m_odeYRows, 1);
2765     pDblY->set(y);
2766     pDblY->IncreaseRef();
2767     in.push_back(pDblY);
2768
2769     types::Double* pDblYdot = new types::Double(m_odeYRows, 1);
2770     pDblYdot->set(ydot);
2771     pDblYdot->IncreaseRef();
2772     in.push_back(pDblYdot);
2773
2774     types::Double* pDblCj = new types::Double(*cj);
2775     pDblCj->IncreaseRef();
2776     in.push_back(pDblCj);
2777
2778     for (int i = 0; i < (int)m_JacArgs.size(); i++)
2779     {
2780         m_JacArgs[i]->IncreaseRef();
2781         in.push_back(m_JacArgs[i]);
2782     }
2783
2784     try
2785     {
2786         // new std::string("") is delete in destructor of ast::CommentExp
2787         m_pCallJacFunction->invoke(in, opt, iRetCount, out, ast::CommentExp(Location(), new std::string("")));
2788     }
2789     catch (const ast::InternalError& ie)
2790     {
2791         for (int i = 0; i < (int)m_JacArgs.size(); i++)
2792         {
2793             m_JacArgs[i]->DecreaseRef();
2794         }
2795
2796         throw ie;
2797     }
2798
2799     for (int i = 0; i < (int)m_JacArgs.size(); i++)
2800     {
2801         m_JacArgs[i]->DecreaseRef();
2802     }
2803
2804     if (out.size() != iRetCount)
2805     {
2806         const char* pstrName = m_pCallJacFunction->getName().c_str();
2807         sprintf(errorMsg, _("%s: Wrong number of output argument(s): %d expected.\n"), pstrName, iRetCount);
2808         throw ast::InternalError(errorMsg);
2809     }
2810
2811     out[0]->IncreaseRef();
2812
2813     pDblT->DecreaseRef();
2814     if (pDblT->isDeletable())
2815     {
2816         delete pDblT;
2817     }
2818
2819     pDblY->DecreaseRef();
2820     if (pDblY->isDeletable())
2821     {
2822         delete pDblY;
2823     }
2824
2825     pDblYdot->DecreaseRef();
2826     if (pDblYdot->isDeletable())
2827     {
2828         delete pDblYdot;
2829     }
2830
2831     pDblCj->DecreaseRef();
2832     if (pDblCj->isDeletable())
2833     {
2834         delete pDblCj;
2835     }
2836
2837     out[0]->DecreaseRef();
2838     if (out[0]->isDouble() == false)
2839     {
2840         const char* pstrName = m_pCallJacFunction->getName().c_str();
2841         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Real matrix expected.\n"), pstrName, 1);
2842         throw ast::InternalError(errorMsg);
2843     }
2844
2845     types::Double* pDblOutPd = out[0]->getAs<types::Double>();
2846     if ( (pDblOutPd->getCols() != m_odeYRows) ||
2847             (!m_bandedJac && pDblOutPd->getRows() != m_odeYRows) ||
2848             (m_bandedJac && pDblOutPd->getRows() != (2 * m_ml + m_mu + 1)))
2849     {
2850         const char* pstrName = m_pCallJacFunction->getName().c_str();
2851         sprintf(errorMsg, _("%s: Wrong size for output argument #%d: A matrix of size %d x %d expected.\n"), pstrName, 1, m_odeYRows, (2 * m_ml + m_mu + 1));
2852         throw ast::InternalError(errorMsg);
2853     }
2854
2855     int size = pDblOutPd->getSize();
2856     C2F(dcopy)(&size, pDblOutPd->get(), &one, pd, &one);
2857     if (out[0]->isDeletable())
2858     {
2859         delete out[0];
2860     }
2861 }
2862
2863 void DifferentialEquationFunctions::callDasrtMacroG(int* ny, double* t, double* y, int* ng, double* gout, double* rpar, int* ipar)
2864 {
2865     char errorMsg[256];
2866     int one         = 1;
2867     int iRetCount   = 1;
2868
2869     types::typed_list in;
2870     types::typed_list out;
2871     types::optional_list opt;
2872
2873     types::Double* pDblT = new types::Double(*t);
2874     pDblT->IncreaseRef();
2875     in.push_back(pDblT);
2876
2877     types::Double* pDblY = new types::Double(*ny, 1);
2878     pDblY->set(y);
2879     pDblY->IncreaseRef();
2880     in.push_back(pDblY);
2881
2882     for (int i = 0; i < (int)m_odeGArgs.size(); i++)
2883     {
2884         m_odeGArgs[i]->IncreaseRef();
2885         in.push_back(m_odeGArgs[i]);
2886     }
2887
2888     try
2889     {
2890         // new std::string("") is delete in destructor of ast::CommentExp
2891         m_pCallGFunction->invoke(in, opt, iRetCount, out, ast::CommentExp(Location(), new std::string("")));
2892     }
2893     catch (const ast::InternalError& ie)
2894     {
2895         for (int i = 0; i < (int)m_odeGArgs.size(); i++)
2896         {
2897             m_odeGArgs[i]->DecreaseRef();
2898         }
2899
2900         throw ie;
2901     }
2902
2903     for (int i = 0; i < (int)m_odeGArgs.size(); i++)
2904     {
2905         m_odeGArgs[i]->DecreaseRef();
2906     }
2907
2908     if (out.size() != iRetCount)
2909     {
2910         const char* pstrName = m_pCallGFunction->getName().c_str();
2911         sprintf(errorMsg, _("%s: Wrong number of output argument(s): %d expected.\n"), pstrName, iRetCount);
2912         throw ast::InternalError(errorMsg);
2913     }
2914
2915     out[0]->IncreaseRef();
2916
2917     pDblT->DecreaseRef();
2918     if (pDblT->isDeletable())
2919     {
2920         delete pDblT;
2921     }
2922
2923     pDblY->DecreaseRef();
2924     if (pDblY->isDeletable())
2925     {
2926         delete pDblY;
2927     }
2928
2929     out[0]->DecreaseRef();
2930     if (out[0]->isDouble() == false)
2931     {
2932         const char* pstrName = m_pCallGFunction->getName().c_str();
2933         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Real matrix expected.\n"), pstrName, 1);
2934         throw ast::InternalError(errorMsg);
2935     }
2936
2937     types::Double* pDblOutGout = out[0]->getAs<types::Double>();
2938     if (pDblOutGout->getSize() != *ng)
2939     {
2940         const char* pstrName = m_pCallGFunction->getName().c_str();
2941         sprintf(errorMsg, _("%s: Wrong size for output argument #%d: A matrix of size %d expected.\n"), pstrName, 1, *ng);
2942         throw ast::InternalError(errorMsg);
2943     }
2944
2945     C2F(dcopy)(ng, pDblOutGout->get(), &one, gout, &one);
2946     if (out[0]->isDeletable())
2947     {
2948         delete out[0];
2949     }
2950 }
2951
2952 void DifferentialEquationFunctions::callDaskrMacroPjac(double* res, int* ires, int* neq, double* t, double* y, double* ydot,
2953         double* rewt, double* savr, double* wk, double* h, double* cj,
2954         double* wp, int* iwp, int* ier, double* rpar, int* ipar)
2955 {
2956     // macro : [R, iR, ier] = psol(neq, t, y, ydot, h, cj, rewt, savr)
2957     char errorMsg[256];
2958     int one         = 1;
2959     int iRetCount   = 3;
2960
2961     types::typed_list in;
2962     types::typed_list out;
2963     types::optional_list opt;
2964
2965     types::Double* pDblNeq = new types::Double((double)(*neq));
2966     pDblNeq->IncreaseRef();
2967     in.push_back(pDblNeq);
2968
2969     types::Double* pDblT = new types::Double(*t);
2970     pDblT->IncreaseRef();
2971     in.push_back(pDblT);
2972
2973     types::Double* pDblY = new types::Double(m_odeYRows, 1);
2974     pDblY->set(y);
2975     pDblY->IncreaseRef();
2976     in.push_back(pDblY);
2977
2978     types::Double* pDblYdot = new types::Double(m_odeYRows, 1);
2979     pDblYdot->set(ydot);
2980     pDblYdot->IncreaseRef();
2981     in.push_back(pDblYdot);
2982
2983     types::Double* pDblH = new types::Double(*h);
2984     pDblH->IncreaseRef();
2985     in.push_back(pDblH);
2986
2987     types::Double* pDblCj = new types::Double(*cj);
2988     pDblCj->IncreaseRef();
2989     in.push_back(pDblCj);
2990
2991     types::Double* pDblRewt = new types::Double(m_odeYRows, 1);
2992     pDblRewt->set(rewt);
2993     pDblRewt->IncreaseRef();
2994     in.push_back(pDblRewt);
2995
2996     types::Double* pDblSavr = new types::Double(m_odeYRows, 1);
2997     pDblSavr->set(savr);
2998     pDblSavr->IncreaseRef();
2999     in.push_back(pDblSavr);
3000
3001     for (int i = 0; i < (int)m_pJacArgs.size(); i++)
3002     {
3003         m_pJacArgs[i]->IncreaseRef();
3004         in.push_back(m_pJacArgs[i]);
3005     }
3006
3007     try
3008     {
3009         // new std::string("") is delete in destructor of ast::CommentExp
3010         m_pCallPjacFunction->invoke(in, opt, iRetCount, out, ast::CommentExp(Location(), new std::string("")));
3011     }
3012     catch (const ast::InternalError& ie)
3013     {
3014         for (int i = 0; i < (int)m_pJacArgs.size(); i++)
3015         {
3016             m_pJacArgs[i]->DecreaseRef();
3017         }
3018
3019         throw ie;
3020     }
3021
3022     for (int i = 0; i < (int)m_pJacArgs.size(); i++)
3023     {
3024         m_pJacArgs[i]->DecreaseRef();
3025     }
3026
3027     if (out.size() != iRetCount)
3028     {
3029         const char* pstrName = m_pCallPjacFunction->getName().c_str();
3030         sprintf(errorMsg, _("%s: Wrong number of output argument(s): %d expected.\n"), pstrName, iRetCount);
3031         throw ast::InternalError(errorMsg);
3032     }
3033
3034     out[0]->IncreaseRef();
3035     out[1]->IncreaseRef();
3036     out[2]->IncreaseRef();
3037
3038     pDblNeq->DecreaseRef();
3039     if (pDblNeq->isDeletable())
3040     {
3041         delete pDblNeq;
3042     }
3043
3044     pDblT->DecreaseRef();
3045     if (pDblT->isDeletable())
3046     {
3047         delete pDblT;
3048     }
3049
3050     pDblY->DecreaseRef();
3051     if (pDblY->isDeletable())
3052     {
3053         delete pDblY;
3054     }
3055
3056     pDblYdot->DecreaseRef();
3057     if (pDblYdot->isDeletable())
3058     {
3059         delete pDblYdot;
3060     }
3061
3062     pDblH->DecreaseRef();
3063     if (pDblH->isDeletable())
3064     {
3065         delete pDblH;
3066     }
3067
3068     pDblCj->DecreaseRef();
3069     if (pDblCj->isDeletable())
3070     {
3071         delete pDblCj;
3072     }
3073
3074     pDblRewt->DecreaseRef();
3075     if (pDblRewt->isDeletable())
3076     {
3077         delete pDblRewt;
3078     }
3079
3080     pDblSavr->DecreaseRef();
3081     if (pDblSavr->isDeletable())
3082     {
3083         delete pDblSavr;
3084     }
3085
3086     // check type of output arguments
3087     if (out[0]->isDouble() == false)
3088     {
3089         const char* pstrName = m_pCallPjacFunction->getName().c_str();
3090         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Real matrix expected.\n"), pstrName, 1);
3091         throw ast::InternalError(errorMsg);
3092     }
3093
3094     if (out[1]->isDouble() == false)
3095     {
3096         const char* pstrName = m_pCallPjacFunction->getName().c_str();
3097         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Real matrix expected.\n"), pstrName, 2);
3098         throw ast::InternalError(errorMsg);
3099     }
3100
3101     if (out[2]->isDouble() == false)
3102     {
3103         const char* pstrName = m_pCallPjacFunction->getName().c_str();
3104         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Real matrix expected.\n"), pstrName, 3);
3105         throw ast::InternalError(errorMsg);
3106     }
3107
3108     //  return [R, iR, ier]
3109     types::Double* pDblOutWp  = out[0]->getAs<types::Double>();
3110     types::Double* pDblOutIwp = out[1]->getAs<types::Double>();
3111     types::Double* pDblOutIer = out[2]->getAs<types::Double>();
3112
3113     // check size of output argument
3114     if (pDblOutWp->getSize() != *neq **neq)
3115     {
3116         const char* pstrName = m_pCallPjacFunction->getName().c_str();
3117         sprintf(errorMsg, _("%s: Wrong size for output argument #%d: A matrix of size %d expected.\n"), pstrName, 1, *neq **neq);
3118         throw ast::InternalError(errorMsg);
3119     }
3120
3121     if (pDblOutIwp->getSize() != 2 * *neq **neq)
3122     {
3123         const char* pstrName = m_pCallPjacFunction->getName().c_str();
3124         sprintf(errorMsg, _("%s: Wrong size for output argument #%d: A matrix of size %d expected.\n"), pstrName, 2, 2 * *neq **neq);
3125         throw ast::InternalError(errorMsg);
3126     }
3127
3128     if (pDblOutIer->isScalar() == false)
3129     {
3130         const char* pstrName = m_pCallPjacFunction->getName().c_str();
3131         sprintf(errorMsg, _("%s: Wrong size for output argument #%d: A Scalar expected.\n"), pstrName, 3);
3132         throw ast::InternalError(errorMsg);
3133     }
3134
3135     // copy output macro results in output variables
3136     int size = pDblOutWp->getSize();
3137     C2F(dcopy)(&size, pDblOutWp->get(), &one, wp, &one);
3138
3139     double* pdblIwp = pDblOutIwp->get();
3140     for (int i = 0; i < pDblOutIwp->getSize(); i++)
3141     {
3142         iwp[i] = (int)pdblIwp[i];
3143     }
3144
3145     *ier = (int)(pDblOutIer->get(0));
3146
3147     // delete output macro result
3148     out[0]->DecreaseRef();
3149     if (out[0]->isDeletable())
3150     {
3151         delete out[0];
3152     }
3153
3154     out[1]->DecreaseRef();
3155     if (out[1]->isDeletable())
3156     {
3157         delete out[1];
3158     }
3159
3160     out[2]->DecreaseRef();
3161     if (out[2]->isDeletable())
3162     {
3163         delete out[2];
3164     }
3165 }
3166
3167 void DifferentialEquationFunctions::callDaskrMacroPsol(int* neq, double* t, double* y, double* ydot, double* savr, double* wk,
3168         double* cj, double* wght, double* wp, int* iwp, double* b, double* eplin,
3169         int* ier, double* rpar, int* ipar)
3170 {
3171     // macro : [b, ier] = psol(R, iR, b)
3172     char errorMsg[256];
3173     int one         = 1;
3174     int iRetCount   = 2;
3175
3176     types::typed_list in;
3177     types::typed_list out;
3178     types::optional_list opt;
3179
3180     // input arguments psol(R, iR, b)
3181     types::Double* pDblR = new types::Double(*neq **neq, 1);
3182     pDblR->set(wp);
3183     pDblR->IncreaseRef();
3184     in.push_back(pDblR);
3185
3186     types::Double* pDblIR = new types::Double(*neq **neq, 2);
3187     double* pdblIR = pDblIR->get();
3188     for (int i = 0; i < pDblIR->getSize(); i++)
3189     {
3190         pdblIR[i] = (double)iwp[i];
3191     }
3192     pDblIR->IncreaseRef();
3193     in.push_back(pDblIR);
3194
3195     types::Double* pDblB = new types::Double(*neq, 1);
3196     pDblB->set(b);
3197     pDblB->IncreaseRef();
3198     in.push_back(pDblB);
3199
3200     // optional arguments
3201     for (int i = 0; i < (int)m_pSolArgs.size(); i++)
3202     {
3203         m_pSolArgs[i]->IncreaseRef();
3204         in.push_back(m_pSolArgs[i]);
3205     }
3206
3207     try
3208     {
3209         // new std::string("") is delete in destructor of ast::CommentExp
3210         m_pCallPsolFunction->invoke(in, opt, iRetCount, out, ast::CommentExp(Location(), new std::string("")));
3211     }
3212     catch (const ast::InternalError& ie)
3213     {
3214         for (int i = 0; i < (int)m_pSolArgs.size(); i++)
3215         {
3216             m_pSolArgs[i]->DecreaseRef();
3217         }
3218
3219         throw ie;
3220     }
3221
3222     for (int i = 0; i < (int)m_pSolArgs.size(); i++)
3223     {
3224         m_pSolArgs[i]->DecreaseRef();
3225     }
3226
3227     // get output
3228     if (out.size() != iRetCount)
3229     {
3230         const char* pstrName = m_pCallPsolFunction->getName().c_str();
3231         sprintf(errorMsg, _("%s: Wrong number of output argument(s): %d expected.\n"), pstrName, iRetCount);
3232         throw ast::InternalError(errorMsg);
3233     }
3234
3235     out[0]->IncreaseRef();
3236     out[1]->IncreaseRef();
3237
3238     // free input arguments
3239     pDblR->DecreaseRef();
3240     if (pDblR->isDeletable())
3241     {
3242         delete pDblR;
3243     }
3244
3245     pDblIR->DecreaseRef();
3246     if (pDblIR->isDeletable())
3247     {
3248         delete pDblIR;
3249     }
3250
3251     pDblB->DecreaseRef();
3252     if (pDblB->isDeletable())
3253     {
3254         delete pDblB;
3255     }
3256
3257     // check output result
3258     if (out[0]->isDouble() == false)
3259     {
3260         const char* pstrName = m_pCallPsolFunction->getName().c_str();
3261         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Real matrix expected.\n"), pstrName, 1);
3262         throw ast::InternalError(errorMsg);
3263     }
3264
3265     if (out[1]->isDouble() == false)
3266     {
3267         const char* pstrName = m_pCallPsolFunction->getName().c_str();
3268         sprintf(errorMsg, _("%s: Wrong type for output argument #%d: Real matrix expected.\n"), pstrName, 2);
3269         throw ast::InternalError(errorMsg);
3270     }
3271
3272     // return arguments [b, ier] = psol()
3273     types::Double* pDblOutB  = out[0]->getAs<types::Double>();
3274     if (pDblOutB->getSize() != *neq) // size of b is neq
3275     {
3276         const char* pstrName = m_pCallPsolFunction->getName().c_str();
3277         sprintf(errorMsg, _("%s: Wrong size for output argument #%d: A matrix of size %d expected.\n"), pstrName, 1, *neq);
3278         throw ast::InternalError(errorMsg);
3279     }
3280
3281     // get scalar ier
3282     types::Double* pDblOutIer = out[1]->getAs<types::Double>();
3283     if (pDblOutIer->isScalar() == false)
3284     {
3285         const char* pstrName = m_pCallPsolFunction->getName().c_str();
3286         sprintf(errorMsg, _("%s: Wrong size for output argument #%d: A Scalar expected.\n"), pstrName, 2);
3287         throw ast::InternalError(errorMsg);
3288     }
3289
3290     // copy output macro results in output variables
3291     C2F(dcopy)(neq, pDblOutB->get(), &one, b, &one);
3292     *ier = (int)(pDblOutIer->get(0));
3293
3294     // free output arguments
3295     out[0]->DecreaseRef();
3296     if (out[0]->isDeletable())
3297     {
3298         delete out[0];
3299     }
3300
3301     out[1]->DecreaseRef();
3302     if (out[1]->isDeletable())
3303     {
3304         delete out[1];
3305     }
3306 }