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