2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2011 - DIGITEO - Cedric DELAMARRE
5 * Copyright (C) 2012 - 2016 - Scilab Enterprises
7 * This file is hereby licensed under the terms of the GNU GPL v2.0,
8 * pursuant to article 5.3.4 of the CeCILL v.2.1.
9 * This file was originally licensed under the terms of the CeCILL v2.1,
10 * and continues to be available under such terms.
11 * For more information, see the COPYING file which you should have received
12 * along with this program.
15 /*--------------------------------------------------------------------------*/
17 #include "differential_equations_gw.hxx"
18 #include "function.hxx"
22 #include "callable.hxx"
23 #include "differentialequationfunctions.hxx"
24 #include "runvisitor.hxx"
25 #include "context.hxx"
26 #include "checkodeerror.hxx"
30 #include "localization.h"
32 #include "scifunctions.h"
33 #include "elem_common.h"
34 #include "configvariable_interface.h"
36 #include "common_structure.h"
37 #include "sci_malloc.h"
40 /*--------------------------------------------------------------------------*/
42 types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, types::typed_list &out)
45 types::String* pStrType = NULL;
46 const wchar_t* wcsType = L"lsoda";
50 // types::Polynom* pPolyY0 = NULL;
51 types::Double* pDblY0 = NULL;
52 double* pdYData = NULL; // contain y0 following by all args data in list case.
53 int sizeOfpdYData = 0;
56 types::Double* pDblT0 = NULL;
57 types::Double* pDblT = NULL;
58 types::Double* pDblRtol = NULL;
59 types::Double* pDblAtol = NULL;
60 types::Double* pDblNg = NULL;
61 types::Double* pDblW = NULL;
62 types::Double* pDblIw = NULL;
65 types::Double* pDblOdeOptions = NULL;
68 types::Double* pDblYOut = NULL;
69 // types::Double* pDblTOut = NULL;
70 // types::Polynom* pPolyYOut = NULL;
72 // Indicate if the function is given.
74 bool bFuncJac = false;
77 int iPos = 0; // Position in types::typed_list in
80 int* YSize = NULL; // YSize(1) = size of y0,
81 // YSize(n) = size of Args(n) in list case.
83 C2F(eh0001).mesflg = 1; // flag to control printing of error messages in lapack routine.
84 // 1 means print, 0 means no printing.
85 C2F(eh0001).lunit = 6; // 6 = stdout
87 int one = 1; // used in C2F(dcopy)
92 // error message catched
93 std::wostringstream os;
96 // *** check the minimal number of input args. ***
99 Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ode", 4);
100 return types::Function::Error;
103 // *** Get the methode. ***
104 if (in[0]->isString())
106 pStrType = in[0]->getAs<types::String>();
107 wcsType = pStrType->get(0);
113 if (wcscmp(wcsType, L"adams") == 0)
117 else if (wcscmp(wcsType, L"stiff") == 0)
121 else if (wcscmp(wcsType, L"root") == 0)
125 else if (wcscmp(wcsType, L"roots") == 0)
127 sciprint(_("%s: Feature %s is obsolete.\n"), _("Warning"), "roots");
128 sciprint(_("%s: Please use %s instead.\n"), _("Warning"), "root");
131 else if (wcscmp(wcsType, L"discrete") == 0)
135 else if (wcscmp(wcsType, L"rk") == 0)
139 else if (wcscmp(wcsType, L"rkf") == 0)
143 else if (wcscmp(wcsType, L"fix") == 0)
149 Scierror(999, _("%s: Wrong value for input argument #%d: It must be one of the following strings: adams, stiff, rk, rkf, fix, root or discrete.\n"), "ode", 1);
150 return types::Function::Error;
154 // *** check number of output args according the methode. ***
159 Scierror(78, _("%s: Wrong number of output argument(s): %d or %d expected.\n"), "ode", 1, 3);
160 return types::Function::Error;
165 if (_iRetCount == 3 || _iRetCount > 4)
167 Scierror(78, _("%s: Wrong number of output argument(s): %d, %d or %d expected.\n"), "ode", 1, 2, 4);
168 return types::Function::Error;
175 Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "ode", 1);
176 return types::Function::Error;
180 // *** check type of input args and get it. ***
182 if (in[iPos]->isDouble())
184 pDblY0 = in[iPos]->getAs<types::Double>();
185 if (pDblY0->isComplex())
187 Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "ode", iPos + 1);
188 return types::Function::Error;
192 else if(in[iPos]->isPoly())
194 pPolyY0 = in[iPos]->getAs<types::Polynom>();
195 if(pPolyY0->isComplex())
197 Scierror(999, _("%s: Wrong type for input argument #%d: A real polynom expected.\n"), "ode", iPos+1);
198 return types::Function::Error;
201 if(pPolyY0->isScalar() == false)
203 Scierror(999, _("%s: Wrong size for input argument #%d: A real scalar polynom expected.\n"), "ode", iPos+1);
204 return types::Function::Error;
207 double* dbl = (double*)MALLOC(pPolyY0->getCoef()->getSize() * sizeof(double));
208 vTransposeRealMatrix( pPolyY0->getCoef()->get(),
209 pPolyY0->getCoef()->getRows(),
210 pPolyY0->getCoef()->getCols(),
213 pDblY0 = new types::Double(pPolyY0->getCoef()->getCols(), pPolyY0->getCoef()->getRows());
219 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ode", iPos + 1);
220 return types::Function::Error;
225 if (in[iPos]->isDouble() == false)
227 Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "ode", iPos + 1);
228 return types::Function::Error;
231 pDblT0 = in[iPos]->getAs<types::Double>();
233 if (pDblT0->isScalar() == false)
235 Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "ode", iPos + 1);
236 return types::Function::Error;
241 if (in[iPos]->isDouble() == false)
243 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ode", iPos + 1);
244 return types::Function::Error;
247 pDblT = in[iPos]->getAs<types::Double>();
250 DifferentialEquationFunctions deFunctionsManager(L"ode");
251 DifferentialEquation::addDifferentialEquationFunctions(&deFunctionsManager);
252 deFunctionsManager.setOdeYRows(pDblY0->getRows());
253 deFunctionsManager.setOdeYCols(pDblY0->getCols());
255 YSize = (int*)MALLOC(sizeOfYSize * sizeof(int));
256 *YSize = pDblY0->getSize();
257 pdYData = (double*)MALLOC(pDblY0->getSize() * sizeof(double));
258 C2F(dcopy)(YSize, pDblY0->get(), &one, pdYData, &one);
264 Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ode", 5);
265 DifferentialEquation::removeDifferentialEquationFunctions();
268 return types::Function::Error;
271 if (in[4]->isCallable() == false && in[4]->isString() == false && in[4]->isList() == false)
273 Scierror(999, _("%s: Wrong type for input argument #%d: A function expected.\n"), "ode", 5);
274 DifferentialEquation::removeDifferentialEquationFunctions();
277 return types::Function::Error;
281 for (iPos++; iPos < (int)in.size(); iPos++)
283 if (in[iPos]->isDouble())
285 if (pDblRtol == NULL && bFuncF == false)
287 pDblRtol = in[iPos]->getAs<types::Double>();
288 if (pDblRtol->getSize() != pDblY0->getSize() && pDblRtol->isScalar() == false)
290 Scierror(267, _("%s: Arg %d and arg %d must have equal dimensions.\n"), "ode", pStrType ? 2 : 1, iPos + 1);
291 DifferentialEquation::removeDifferentialEquationFunctions();
294 return types::Function::Error;
297 else if (pDblAtol == NULL && bFuncF == false)
299 pDblAtol = in[iPos]->getAs<types::Double>();
300 if (pDblAtol->getSize() != pDblY0->getSize() && pDblAtol->isScalar() == false)
302 Scierror(267, _("%s: Arg %d and arg %d must have equal dimensions.\n"), "ode", pStrType ? 2 : 1, iPos + 1);
303 DifferentialEquation::removeDifferentialEquationFunctions();
306 return types::Function::Error;
309 else if (pDblNg == NULL && bFuncF == true && meth == 3)
311 pDblNg = in[iPos]->getAs<types::Double>();
313 else if (pDblW == NULL && bFuncF == true && (bFuncG == true || meth != 3))
315 if ((int)in.size() == iPos + 2)
317 if (in[iPos + 1]->isDouble() == false)
319 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ode", iPos + 2);
320 DifferentialEquation::removeDifferentialEquationFunctions();
323 return types::Function::Error;
326 pDblW = in[iPos]->getAs<types::Double>();
327 pDblIw = in[iPos + 1]->getAs<types::Double>();
332 Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ode", iPos + 2);
333 DifferentialEquation::removeDifferentialEquationFunctions();
336 return types::Function::Error;
341 Scierror(999, _("%s: Wrong type for input argument #%d: A function expected.\n"), "ode", iPos + 1);
342 DifferentialEquation::removeDifferentialEquationFunctions();
345 return types::Function::Error;
348 else if (in[iPos]->isCallable())
350 types::Callable* pCall = in[iPos]->getAs<types::Callable>();
353 deFunctionsManager.setFFunction(pCall);
356 else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
358 deFunctionsManager.setJacFunction(pCall);
361 else if (bFuncG == false && meth == 3)
363 deFunctionsManager.setGFunction(pCall);
368 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ode", iPos + 1);
369 DifferentialEquation::removeDifferentialEquationFunctions();
372 return types::Function::Error;
375 else if (in[iPos]->isString())
377 types::String* pStr = in[iPos]->getAs<types::String>();
382 bOK = deFunctionsManager.setFFunction(pStr);
385 else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
387 bOK = deFunctionsManager.setJacFunction(pStr);
390 else if (bFuncG == false && meth == 3)
392 bOK = deFunctionsManager.setGFunction(pStr);
397 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ode", iPos + 1);
398 DifferentialEquation::removeDifferentialEquationFunctions();
401 return types::Function::Error;
406 char* pst = wide_string_to_UTF8(pStr->get(0));
407 Scierror(50, _("%s: Subroutine not found: %s\n"), "ode", pst);
409 DifferentialEquation::removeDifferentialEquationFunctions();
412 return types::Function::Error;
415 else if (in[iPos]->isList())
417 types::List* pList = in[iPos]->getAs<types::List>();
419 if (pList->getSize() == 0)
421 Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "ode", iPos + 1, "(string empty)");
422 DifferentialEquation::removeDifferentialEquationFunctions();
425 return types::Function::Error;
428 if (bFuncF && (bFuncJac || pDblNg) && (bFuncG || meth != 3))
430 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ode", iPos + 1);
431 DifferentialEquation::removeDifferentialEquationFunctions();
434 return types::Function::Error;
437 if (pList->get(0)->isString())
439 types::String* pStr = pList->get(0)->getAs<types::String>();
445 bOK = deFunctionsManager.setFFunction(pStr);
446 sizeOfpdYData = *YSize;
448 else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
451 bOK = deFunctionsManager.setJacFunction(pStr);
452 if (sizeOfpdYData == 0)
454 sizeOfpdYData = *YSize;
457 else if (bFuncG == false && meth == 3)
460 bOK = deFunctionsManager.setGFunction(pStr);
461 if (sizeOfpdYData == 0)
463 sizeOfpdYData = *YSize;
469 char* pst = wide_string_to_UTF8(pStr->get(0));
470 Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "ode", iPos + 1, pst);
472 DifferentialEquation::removeDifferentialEquationFunctions();
475 return types::Function::Error;
478 int* sizeTemp = YSize;
479 int totalSize = sizeOfpdYData;
481 YSize = (int*)MALLOC((sizeOfYSize + pList->getSize() - 1) * sizeof(int));
482 memcpy(YSize, sizeTemp, sizeOfYSize * sizeof(int));
484 std::vector<types::Double*> vpDbl;
485 for (int iter = 0; iter < pList->getSize() - 1; iter++)
487 if (pList->get(iter + 1)->isDouble() == false)
489 Scierror(999, _("%s: Wrong type for input argument #%d: Argument %d in the list must be a matrix.\n"), "ode", iPos + 1, iter + 1);
490 DifferentialEquation::removeDifferentialEquationFunctions();
493 return types::Function::Error;
496 vpDbl.push_back(pList->get(iter + 1)->getAs<types::Double>());
497 YSize[sizeOfYSize + iter] = vpDbl[iter]->getSize();
498 totalSize += YSize[sizeOfYSize + iter];
501 double* pdYDataTemp = pdYData;
502 pdYData = (double*)MALLOC(totalSize * sizeof(double));
503 C2F(dcopy)(&sizeOfpdYData, pdYDataTemp, &one, pdYData, &one);
505 int position = sizeOfpdYData;
506 for (int iter = 0; iter < pList->getSize() - 1; iter++)
508 C2F(dcopy)(&YSize[sizeOfYSize + iter], vpDbl[iter]->get(), &one, &pdYData[position], &one);
509 position += vpDbl[iter]->getSize();
512 sizeOfpdYData = totalSize;
513 sizeOfYSize += pList->getSize() - 1;
517 else if (pList->get(0)->isCallable())
522 deFunctionsManager.setFFunction(pList->get(0)->getAs<types::Callable>());
523 for (int iter = 1; iter < pList->getSize(); iter++)
525 deFunctionsManager.setFArgs(pList->get(iter)->getAs<types::InternalType>());
528 else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
531 deFunctionsManager.setJacFunction(pList->get(0)->getAs<types::Callable>());
532 for (int iter = 1; iter < pList->getSize(); iter++)
534 deFunctionsManager.setJacArgs(pList->get(iter)->getAs<types::InternalType>());
537 else if (bFuncG == false && meth == 3)
540 deFunctionsManager.setGFunction(pList->get(0)->getAs<types::Callable>());
541 for (int iter = 1; iter < pList->getSize(); iter++)
543 deFunctionsManager.setGArgs(pList->get(iter)->getAs<types::InternalType>());
549 Scierror(999, _("%s: Wrong type for input argument #%d: The first argument in the list must be a string or a function.\n"), "ode", iPos + 1);
550 DifferentialEquation::removeDifferentialEquationFunctions();
553 return types::Function::Error;
558 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a function expected.\n"), "ode", iPos + 1);
559 DifferentialEquation::removeDifferentialEquationFunctions();
562 return types::Function::Error;
568 int val = (meth == 3) ? 3 : 1;
569 Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ode", in.size() + val);
570 DifferentialEquation::removeDifferentialEquationFunctions();
573 return types::Function::Error;
575 if (pDblNg == NULL && meth == 3)
577 Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ode", in.size() + 2);
578 DifferentialEquation::removeDifferentialEquationFunctions();
581 return types::Function::Error;
583 if (bFuncG == false && meth == 3)
585 Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ode", in.size() + 1);
586 DifferentialEquation::removeDifferentialEquationFunctions();
589 return types::Function::Error;
592 // *** Initialization. ***
595 int itask = 1; // itask = 1 for normal computation of output values of y at t = tout.
596 int istate = 1; // istate = integer flag (input and output). set istate = 1.
597 int iopt = 0; // iopt = 0 to indicate no optional inputs used.
598 int jt = 2; // jacobian type indicator. set jt = 2.
603 double* rwork = NULL;
608 // contain ls0001, lsa001 and eh0001 structures
609 double* dStructTab = NULL;
610 int* iStructTab = NULL;
611 int dStructTabSize = 0;
612 int iStructTabSize = 0;
614 int rwSize = 0; // rwSize = dStructTab + rworkSize
615 int iwSize = 0; // iwSize = iStructTab + iworkSize
617 // structures used by lsoda and lsode
618 double* ls0001d = &(C2F(ls0001).tret);
619 int* ls0001i = &(C2F(ls0001).illin);
620 double* lsa001d = &(C2F(lsa001).tsw);
621 int* lsa001i = &(C2F(lsa001).insufr);
622 double* lsr001d = &(C2F(lsr001).rownr3[0]);
623 int* lsr001i = &(C2F(lsr001).lg0);
624 int* eh0001i = &(C2F(eh0001).mesflg);
627 types::InternalType* pIT = symbol::Context::getInstance()->get(symbol::Symbol(L"%ODEOPTIONS"));
628 if (pIT != NULL && pIT->isDouble())
630 pDblOdeOptions = pIT->getAs<types::Double>();
631 if (pDblOdeOptions->getSize() == 12)
634 itask = (int)pDblOdeOptions->get(0);
635 jt = (int)pDblOdeOptions->get(5);
636 ml = (int)pDblOdeOptions->get(10);
637 mu = (int)pDblOdeOptions->get(11);
641 if (iopt == 1 && (pDblOdeOptions->get(4) > pDblOdeOptions->get(3))) // hmin > hmax ?
643 Scierror(9999, _("%s: Wrong value of hmin and hmax: hmin = %d is greater than hmax = %d.\n"), "ode", pDblOdeOptions->get(4), pDblOdeOptions->get(3));
644 DifferentialEquation::removeDifferentialEquationFunctions();
647 return types::Function::Error;
650 if (jt < 0 || jt > 5)
652 Scierror(9999, _("%s: Wrong value of Jacobian type: A number between %d and %d expected.\n"), "ode", 0, 5);
653 DifferentialEquation::removeDifferentialEquationFunctions();
656 return types::Function::Error;
659 if (iopt == 0 && bFuncJac)
664 if (bFuncJac && (jt == 2 || jt == 5) && getWarningMode())
666 sciprint(_("%s: Warning: Jacobian is given, but not used.\n"), "ode");
669 if (bFuncJac == false && (jt == 1 || jt == 4))
671 if (getWarningMode())
673 sciprint(_("%s: Warning: No Jacobian external given, but one is required by %ODEOPTIONS(6) value. Jacobian will be estimated.\n"), "ode");
679 //compute itol and set the tolerances rtol and atol.
685 if (pDblRtol->isScalar())
687 rtol = (double*)MALLOC(sizeof(double));
688 *rtol = pDblRtol->get(0);
692 rtol = pDblRtol->get();
698 rtol = (double*)MALLOC(sizeof(double));
699 if (meth == 6 || meth == 7)
711 if (pDblAtol->isScalar())
713 atol = (double*)MALLOC(sizeof(double));
714 *atol = pDblAtol->get(0);
718 atol = pDblAtol->get();
724 atol = (double*)MALLOC(sizeof(double));
725 if (meth == 6 || meth == 7)
735 // Compute rwork, iwork size.
738 if (pDblW) // structure ls0001 have been restored.
740 nyh = C2F(ls0001).nyh;
748 mxordn = (int)pDblOdeOptions->get(7);
749 mxords = (int)pDblOdeOptions->get(8);
752 if (mxordn > 12 || mxords > 5 || mxordn < 1 || mxords < 1)
754 if (getWarningMode())
756 sciprint(_("%s: Warning: Wrong value for maximum stiff/non-stiff order allowed :\nAt most %d for mxordn, %d for mxords and no null value for both expected.\nWrong value will be reduced to the default value.\n"), "ode", 12, 5);
770 case 3: // lsodar (root)
772 // lrn = 20 + nyh*(mxordn+1) + 3*neq + 3*ng,
773 // lrs = 20 + nyh*(mxords+1) + 3*neq + lmat + 3*ng,
775 // nyh = the initial value of neq,
776 // mxordn = 12, unless a smaller value is given as an
778 // mxords = 5, unless a smaller value is given as an
780 // lmat = length of matrix work space..
781 // lmat = neq**2 + 2 if jt = 1 or 2,
782 // lmat = (2*ml + mu + 1)*neq + 2 if jt = 4 or 5.
784 lrn = 3 * (int)pDblNg->get(0);
786 dStructTabSize = 246 - 241;
787 iStructTabSize = 59 - 50;
788 jroot = (int*)MALLOC((int)pDblNg->get(0) * sizeof(int));
792 // lrn = 20 + nyh*(mxordn+1) + 3*neq,
793 // lrs = 20 + nyh*(mxords+1) + 3*neq + lmat,
795 // nyh = the initial value of neq,
796 // mxordn = 12, unless a smaller value is given as an
798 // mxords = 5, unless a smaller value is given as an
800 // lmat = length of matrix work space..
801 // lmat = neq**2 + 2 if jt = 1 or 2,
802 // lmat = (2*ml + mu + 1)*neq + 2 if jt = 4 or 5.
804 if (jt == 1 || jt == 2)
806 lmat = (*YSize) * (*YSize) + 2; // if jt = 1 or 2,
808 else if (jt == 4 || jt == 5)
810 //ml and mu = -1 in all cases
811 lmat = (2 * ml + mu + 1) * (*YSize) + 2; // if jt = 4 or 5.
814 lrn += 20 + nyh * (mxordn + 1) + 3 * (*YSize);
815 lrs += 20 + nyh * (mxords + 1) + 3 * (*YSize) + lmat;
817 rworkSize = std::max(lrn, lrs);
818 iworkSize = 20 + *YSize;
820 dStructTabSize += 241;
821 iStructTabSize += 50;
825 case 1: // lsode (adams) (non stiff)
829 case 2: // lsode (stiff)
831 // 20 + nyh*(maxord + 1) + 3*neq + lmat
833 // nyh = the initial value of neq,
834 // maxord = 12 (if meth = 1) or 5 (if meth = 2) (unless a
835 // smaller value is given as an optional input),
836 // lmat = 0 if miter = 0,
837 // lmat = neq**2 + 2 if miter is 1 or 2,
838 // lmat = neq + 2 if miter = 3, and
839 // lmat = (2*ml+mu+1)*neq + 2 if miter is 4 or 5.
841 if (jt == 1 || jt == 2)
843 lmat = (*YSize) * (*YSize) + 2;
849 else if (jt == 4 || jt == 5)
851 lmat = (2 * ml + mu + 1) * (*YSize) + 2;
854 rworkSize = 20 + nyh * (maxord + 1) + 3 * (*YSize) + lmat;
855 iworkSize = 20; // if jt = 0 or 3
856 if (jt == 1 || jt == 2 || jt == 4 || jt == 5) // iSize = 20 + neq otherwise
858 iworkSize += (*YSize);
861 dStructTabSize = 219;
865 case 4: // lsdisc (discrete)
871 case 5: // lsrgk (rk)
873 rworkSize = 1 + 9 * (*YSize);
877 case 6: // rkf45 (rkf)
879 rworkSize = 3 + 8 * (*YSize);
883 case 7: // rksimp (fix)
885 rworkSize = 3 + 8 * (*YSize);
891 rwSize = rworkSize + dStructTabSize;
892 iwSize = iworkSize + iStructTabSize;
894 rwork = (double*)MALLOC(rworkSize * sizeof(double));
895 iwork = (int*)MALLOC(iworkSize * sizeof(int));
901 if (pDblW->getSize() != rwSize || pDblIw->getSize() != iwSize)
903 Scierror(9999, _("%s: Wrong size for w and iw: w = %d and iw = %d expected.\n"), "ode", rwSize, iwSize);
904 DifferentialEquation::removeDifferentialEquationFunctions();
913 if (itol == 1 || itol == 3)
921 return types::Function::Error;
924 istate = 2; // 1 means this is the first call | 2 means this is not the first call
926 dStructTab = (double*)MALLOC(dStructTabSize * sizeof(double));
927 iStructTab = (int*)MALLOC(iStructTabSize * sizeof(int));
929 C2F(dcopy)(&rworkSize, pDblW->get(), &one, rwork, &one);
930 C2F(dcopy)(&dStructTabSize, pDblW->get() + rworkSize, &one, dStructTab, &one);
932 for (int i = 0; i < iworkSize; i++)
934 iwork[i] = (int)pDblIw->get(i);
937 for (int i = 0; i < iStructTabSize; i++)
939 iStructTab[i] = (int)pDblIw->get(i + iworkSize);
944 //if iopt = 0 default value are used without set it.
949 rwork[0] = pDblOdeOptions->get(1); // tcrit
952 rwork[4] = pDblOdeOptions->get(2); // h0 =0
953 rwork[5] = pDblOdeOptions->get(3); // hmax = %inf
954 rwork[6] = pDblOdeOptions->get(4); // hmin = 0
956 if (jt == 4 || jt == 5)
958 iwork[0] = (int)pDblOdeOptions->get(10);// ml = -1
959 iwork[1] = (int)pDblOdeOptions->get(11);// mu = -1
962 if (meth == 0 || meth == 3)
965 iwork[4] = (int)pDblOdeOptions->get(9); // ixpr = 0
966 iwork[7] = (int)pDblOdeOptions->get(7); // mxordn = 12
967 iwork[8] = (int)pDblOdeOptions->get(8); // mxords = 5
974 iwork[4] = (int)pDblOdeOptions->get(7); // mxordn = 12
978 iwork[4] = (int)pDblOdeOptions->get(8); // mxords = 5
982 iwork[5] = (int)pDblOdeOptions->get(6); // mxstep = 500
983 iwork[6] = 10; // mxhnil = 10 maximum number of messages printed per problem
989 if (meth == 5 || meth == 6)
991 int iLastT = pDblT->getSize() - 1;
992 double t = pDblT->get(iLastT);
993 double t0 = pDblT0->get(0);
994 rwork[(*YSize)] = (t - t0) / 100; // h0
998 if (pDblOdeOptions && pDblOdeOptions->get(9) == 1)
1000 sciprint(_("itask = %d\tmeth = %d\tjactyp = %d\tml = %d\tmu = %d\tiopt = %d\n"), itask, meth, jt, ml, mu, iopt);
1001 sciprint(_("tcrit = %lf\th0 = %lf\thmax = %lf\thmin = %lf\n"), pDblOdeOptions->get(1), pDblOdeOptions->get(2), pDblOdeOptions->get(3), pDblOdeOptions->get(4));
1006 if (pDblW && pDblIw)
1012 //restore ls0001 and eh0001 from w (rwork) and iw (iwork).
1013 C2F(dcopy)(&dSize, dStructTab, &one, ls0001d, &one);
1014 memcpy(ls0001i, iStructTab, 39 * sizeof(int));
1019 //restore lsa001 from w (rwork) and iw (iwork).
1020 if (meth == 0 || meth == 3)
1023 C2F(dcopy)(&dSize, dStructTab + dPos, &one, lsa001d, &one);
1024 memcpy(lsa001i, &iStructTab[iPos], 9 * sizeof(int));
1034 C2F(dcopy)(&dSize, dStructTab + dPos, &one, lsr001d, &one);
1035 memcpy(lsr001i, &iStructTab[iPos], 9 * sizeof(int));
1040 // *** Perform operation. ***
1042 double t0 = pDblT0->get(0);
1043 bool bOneStep = false;
1045 if (itask == 5 || itask == 3 || itask == 2)
1048 if (getWarningMode() && pDblT->isScalar() == false)
1050 sciprint(_("itask = %d: At most one value of t is allowed, the last element of t is used.\n"), itask);
1056 std::list<std::vector<double>> pDblYOutList = std::list<std::vector<double>>();
1057 std::list<double> pDblTOutList = std::list<double>();
1058 int iLastT = pDblT->getSize() - 1;
1059 double t = pDblT->get(iLastT);
1060 int iDir = t - t0 < 0 ? -1 : 1;
1064 std::string strMeth;
1069 case 1: // lsode (adams)
1070 case 2: // lsode (stiff)
1073 int jacType = 10 * meth + jt;
1074 C2F(lsode)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jacType);
1080 int ng = (int)pDblNg->get(0);
1081 C2F(lsodar)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt, ode_g, &ng, jroot);
1084 case 4: // lsdisc (discrete)
1087 C2F(lsdisc)(ode_f, YSize, pdYData, &t0, &t, rwork, &rworkSize, &istate);
1090 case 5: // lsrgk (rk)
1093 C2F(lsrgk)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1096 case 6: // rkf45 (rkf)
1099 C2F(rkf45)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1102 case 7: // rksimp (fix)
1105 C2F(rksimp)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1108 default: // case 0: // lsoda
1111 C2F(lsoda)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt);
1116 err = checkOdeError(meth, istate);
1119 Scierror(999, _("%s: %s exit with state %d.\n"), "ode", strMeth.c_str(), istate);
1122 catch (ast::InternalError &ie)
1124 os << ie.GetErrorMessage();
1129 // FREE allocated data
1130 if (err == 1) // error case
1132 DifferentialEquation::removeDifferentialEquationFunctions();
1149 if (itol == 1 || itol == 3)
1160 wchar_t szError[bsiz];
1161 wchar_t* tmp = to_wide_string(strMeth.c_str());
1162 os_swprintf(szError, bsiz, _W("%ls: An error occurred in '%ls' subroutine.\n").c_str(), L"ode", tmp);
1165 throw ast::InternalError(os.str());
1168 return types::Function::Error;
1171 std::vector<double> newY (pdYData, pdYData + pDblY0->getSize());
1172 pDblYOutList.push_back(newY);
1173 pDblTOutList.push_back(t0);
1175 if (err == 2) // warning case
1177 if (getWarningMode())
1179 sciprint(_("Integration was stopped at t = %lf.\n"), t0);
1184 if (meth == 3 && istate == 3)
1186 // istate == 3 means the integration was successful, and one or more
1187 // roots were found before satisfying the stop condition
1188 // specified by itask. see jroot.
1189 if (getWarningMode())
1191 sciprint(_("%s: Warning: At t = %lf, y is a root, jroot = "), "ode", t0);
1192 for (int k = 0; k < pDblNg->get(0); k++)
1194 sciprint("\t%d", jroot[k]);
1202 while ((t0 - t) * iDir < 0);
1204 int iSizeList = (int)pDblYOutList.size();
1206 if(pPolyY0) // pPolyY0 is scalar
1208 // result format = [2 x iSizeList]
1209 // first lime is the time of result stored in second line.
1210 int size = 2 * iSizeList;
1211 int* iRanks = (int*)MALLOC(size * sizeof(int));
1212 int iMaxRank = pPolyY0->getMaxRank();
1213 for(int i = 0; i < iSizeList; i++)
1215 iRanks[i * 2] = 1; // time rank
1216 iRanks[i * 2 + 1] = iMaxRank; // result rank
1219 pPolyYOut = new types::Polynom(pPolyY0->getVariableName(), 2, iSizeList, iRanks);
1220 pDblYOut = new types::Double(pDblY0->getRows(), 1);
1221 pDblTOut = new types::Double(1, 1);
1223 for(int i = 0; i < iSizeList; i++)
1225 // pDblY0 contain pPolyY0 (ie: pPolyY0 = 2 + x => pDblY0 = [2 1])
1226 for(int j = 0; j < pDblY0->getRows(); j++)
1228 pDblYOut->set(j, pDblYOutList.front()[j]);
1230 pDblTOut->set(0, pDblTOutList.front());
1232 pPolyYOut->setCoef(i * 2, pDblTOut);
1233 pPolyYOut->setCoef(i * 2 + 1, pDblYOut);
1235 pDblYOutList.pop_front();
1236 pDblTOutList.pop_front();
1242 pDblYOut = new types::Double(pDblY0->getRows() + 1, iSizeList);
1243 for (int i = 0; i < iSizeList; i++)
1245 pDblYOut->set(i * (*YSize + 1), pDblTOutList.front());
1246 for (int j = 0; j < *YSize; j++)
1248 pDblYOut->set(i * (*YSize + 1) + (j + 1), pDblYOutList.front()[j]);
1250 pDblYOutList.pop_front();
1251 pDblTOutList.pop_front();
1260 int size = pDblT->getSize();
1261 int* iRanks = (int*)MALLOC(size * sizeof(int));
1262 int iMaxRank = pPolyY0->getMaxRank();
1263 for(int i = 0; i < size; i++)
1265 iRanks[i] = iMaxRank;
1268 pPolyYOut = new types::Polynom(pPolyY0->getVariableName(), 1, pDblT->getSize(), iRanks);
1270 pDblYOut = new types::Double(pDblY0->getRows(), 1);
1275 pDblYOut = new types::Double(pDblY0->getRows(), pDblY0->getCols() * pDblT->getSize());
1277 bool bBreak = false;
1278 for (int i = 0; i < pDblT->getSize(); i++)
1280 double t = pDblT->get(i);
1281 std::string strMeth;
1283 if (itask >= 4 && t > rwork[0]) // rwork[0] => tcrit
1293 case 1: // lsode (adams)
1294 case 2: // lsode (stiff)
1297 int jacType = 10 * meth + jt;
1298 C2F(lsode)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jacType);
1304 int ng = (int)pDblNg->get(0);
1305 C2F(lsodar)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt, ode_g, &ng, jroot);
1308 case 4: // lsdisc (discrete)
1311 C2F(lsdisc)(ode_f, YSize, pdYData, &t0, &t, rwork, &rworkSize, &istate);
1314 case 5: // lsrgk (rk)
1317 C2F(lsrgk)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1320 case 6: // rkf45 (rkf)
1323 C2F(rkf45)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1326 case 7: // rksimp (fix)
1329 C2F(rksimp)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1332 default: // case 0: // lsoda
1335 C2F(lsoda)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt);
1340 err = checkOdeError(meth, istate);
1343 Scierror(999, _("%s: %s exit with state %d.\n"), "ode", strMeth.c_str(), istate);
1346 catch (ast::InternalError &ie)
1348 os << ie.GetErrorMessage();
1352 // FREE allocated data
1353 if (err == 1) // error case
1355 DifferentialEquation::removeDifferentialEquationFunctions();
1372 if (itol == 1 || itol == 3)
1385 wchar_t szError[bsiz];
1386 wchar_t* tmp = to_wide_string(strMeth.c_str());
1387 os_swprintf(szError, bsiz, _W("%ls: An error occurred in '%ls' subroutine.\n").c_str(), L"ode", tmp);
1390 throw ast::InternalError(os.str());
1393 return types::Function::Error;
1398 for(int j = 0; j < pDblY0->getRows(); j++)
1400 pDblYOut->set(j, pdYData[j]);
1402 pPolyYOut->setCoef(i, pDblYOut);
1407 for (int j = 0; j < *YSize; j++)
1409 pDblYOut->set(i * (*YSize) + j, pdYData[j]);
1413 if (err == 2) // warning case
1415 if (getWarningMode())
1417 sciprint(_("Integration was stopped at t = %lf.\n"), t0);
1420 types::Double* pDblYOutTemp = pDblYOut;
1421 pDblYOut = new types::Double(pDblYOutTemp->getRows(), i + 1);
1422 for (int k = 0; k < pDblYOut->getSize(); k++)
1424 pDblYOut->set(k, pDblYOutTemp->get(k));
1426 delete pDblYOutTemp;
1431 if (meth == 3 && istate == 3)
1433 // istate == 3 means the integration was successful, and one or more
1434 // roots were found before satisfying the stop condition
1435 // specified by itask. see jroot.
1437 if (getWarningMode())
1439 sciprint(_("%s: Warning: At t = %lf, y is a root, jroot = "), "ode", t0);
1440 for (int k = 0; k < pDblNg->get(0); k++)
1442 sciprint("\t%d", jroot[k]);
1447 types::Double* pDblYOutTemp = pDblYOut;
1448 pDblYOut = new types::Double(pDblYOutTemp->getRows(), i + 1);
1449 for (int k = 0; k < pDblYOut->getSize(); k++)
1451 pDblYOut->set(k, pDblYOutTemp->get(k));
1453 delete pDblYOutTemp;
1457 if (itask >= 4 && bBreak)
1459 types::Double* pDblYOutTemp = pDblYOut;
1460 pDblYOut = new types::Double(pDblYOutTemp->getRows(), i + 1);
1461 for (int k = 0; k < pDblYOut->getSize(); k++)
1463 pDblYOut->set(k, pDblYOutTemp->get(k));
1465 delete pDblYOutTemp;
1473 if (_iRetCount > 2) //save ls0001 and eh0001 following w and iw.
1479 if (dStructTab == NULL)
1481 dStructTab = (double*)MALLOC(dStructTabSize * sizeof(double));
1482 iStructTab = (int*)MALLOC(iStructTabSize * sizeof(int));
1485 C2F(dcopy)(&dSize, ls0001d, &one, dStructTab, &one);
1486 memcpy(iStructTab, ls0001i, 39 * sizeof(int));
1492 if (meth == 0 || meth == 3)
1495 C2F(dcopy)(&dSize, lsa001d, &one, dStructTab + dPos, &one);
1496 memcpy(&iStructTab[iPos], lsa001i, 9 * sizeof(int));
1506 C2F(dcopy)(&dSize, lsr001d, &one, dStructTab + dPos, &one);
1507 memcpy(&iStructTab[iPos], lsr001i, 9 * sizeof(int));
1512 memcpy(&iStructTab[iPos], eh0001i, 2 * sizeof(int));
1516 // *** Return result in Scilab. ***
1520 out.push_back(pPolyYOut); // y
1525 out.push_back(pDblYOut); // y
1528 if (meth == 3 && _iRetCount >= 2)
1533 for (int i = 0; i < pDblNg->get(0); i++)
1541 types::Double* pDblRd = NULL;
1542 if (sizeOfRd == 1) // Not root found, return empty matrix
1544 pDblRd = types::Double::Empty();
1548 pDblRd = new types::Double(1, sizeOfRd);
1549 //rd: The first entry contains the stopping time.
1550 pDblRd->set(0, C2F(lsr001).tlast);
1551 for (int i = 0; i < pDblNg->get(0); i++)
1556 pDblRd->set(k, (double)i + 1);
1560 out.push_back(pDblRd); // rd
1563 if ((meth < 3 && _iRetCount == 3) || (meth == 3 && _iRetCount == 4))
1565 types::Double* pDblWOut = new types::Double(1, rwSize);
1566 C2F(dcopy)(&rworkSize, rwork, &one, pDblWOut->get(), &one);
1567 C2F(dcopy)(&dStructTabSize, dStructTab, &one, pDblWOut->get() + rworkSize, &one);
1569 types::Double* pDblIwOut = new types::Double(1, iwSize);
1570 for (int i = 0; i < iworkSize; i++)
1572 pDblIwOut->set(i, (double)iwork[i]);
1575 for (int i = 0; i < iStructTabSize; i++)
1577 pDblIwOut->set(iworkSize + i, (double)iStructTab[i]);
1580 out.push_back(pDblWOut); // w
1581 out.push_back(pDblIwOut); // iw
1585 if (itol == 1 || itol == 3) // atol is scalar
1590 if (itol < 3) // rtol is scalar
1611 DifferentialEquation::removeDifferentialEquationFunctions();
1613 return types::Function::OK;