2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2011 - DIGITEO - Cedric DELAMARRE
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
12 /*--------------------------------------------------------------------------*/
14 #include "differential_equations_gw.hxx"
15 #include "function.hxx"
19 #include "callable.hxx"
20 #include "differentialequationfunctions.hxx"
21 #include "runvisitor.hxx"
22 #include "context.hxx"
23 #include "checkodeerror.hxx"
27 #include "localization.h"
29 #include "scifunctions.h"
30 #include "elem_common.h"
31 #include "configvariable_interface.h"
33 #include "common_structure.h"
34 #include "sci_malloc.h"
36 /*--------------------------------------------------------------------------*/
38 types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, types::typed_list &out)
41 types::String* pStrType = NULL;
42 const wchar_t* wcsType = L"lsoda";
46 // types::Polynom* pPolyY0 = NULL;
47 types::Double* pDblY0 = NULL;
48 double* pdYData = NULL; // contain y0 following by all args data in list case.
49 int sizeOfpdYData = 0;
52 types::Double* pDblT0 = NULL;
53 types::Double* pDblT = NULL;
54 types::Double* pDblRtol = NULL;
55 types::Double* pDblAtol = NULL;
56 types::Double* pDblNg = NULL;
57 types::Double* pDblW = NULL;
58 types::Double* pDblIw = NULL;
61 types::Double* pDblOdeOptions = NULL;
64 types::Double* pDblYOut = NULL;
65 // types::Double* pDblTOut = NULL;
66 // types::Polynom* pPolyYOut = NULL;
68 // Indicate if the function is given.
70 bool bFuncJac = false;
73 int iPos = 0; // Position in types::typed_list in
76 int* YSize = NULL; // YSize(1) = size of y0,
77 // YSize(n) = size of Args(n) in list case.
79 C2F(eh0001).mesflg = 1; // flag to control printing of error messages in lapack routine.
80 // 1 means print, 0 means no printing.
81 C2F(eh0001).lunit = 6; // 6 = stdout
83 int one = 1; // used in C2F(dcopy)
88 // *** check the minimal number of input args. ***
91 Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ode", 4);
92 return types::Function::Error;
95 // *** Get the methode. ***
96 if (in[0]->isString())
98 pStrType = in[0]->getAs<types::String>();
99 wcsType = pStrType->get(0);
105 if (wcscmp(wcsType, L"adams") == 0)
109 else if (wcscmp(wcsType, L"stiff") == 0)
113 else if (wcscmp(wcsType, L"root") == 0)
117 else if (wcscmp(wcsType, L"discrete") == 0)
121 else if (wcscmp(wcsType, L"rk") == 0)
125 else if (wcscmp(wcsType, L"rkf") == 0)
129 else if (wcscmp(wcsType, L"fix") == 0)
135 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);
136 return types::Function::Error;
140 // *** check number of output args according the methode. ***
143 if (_iRetCount != 1 && _iRetCount != 3)
145 Scierror(78, _("%s: Wrong number of output argument(s): %d or %d expected.\n"), "ode", 1, 3);
146 return types::Function::Error;
151 if (_iRetCount == 3 || _iRetCount > 4)
153 Scierror(78, _("%s: Wrong number of output argument(s): %d, %d or %d expected.\n"), "ode", 1, 2, 4);
154 return types::Function::Error;
161 Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "ode", 1);
162 return types::Function::Error;
166 // *** check type of input args and get it. ***
168 if (in[iPos]->isDouble())
170 pDblY0 = in[iPos]->getAs<types::Double>();
171 if (pDblY0->isComplex())
173 Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "ode", iPos + 1);
174 return types::Function::Error;
178 else if(in[iPos]->isPoly())
180 pPolyY0 = in[iPos]->getAs<types::Polynom>();
181 if(pPolyY0->isComplex())
183 Scierror(999, _("%s: Wrong type for input argument #%d: A real polynom expected.\n"), "ode", iPos+1);
184 return types::Function::Error;
187 if(pPolyY0->isScalar() == false)
189 Scierror(999, _("%s: Wrong size for input argument #%d: A real scalar polynom expected.\n"), "ode", iPos+1);
190 return types::Function::Error;
193 double* dbl = (double*)MALLOC(pPolyY0->getCoef()->getSize() * sizeof(double));
194 vTransposeRealMatrix( pPolyY0->getCoef()->get(),
195 pPolyY0->getCoef()->getRows(),
196 pPolyY0->getCoef()->getCols(),
199 pDblY0 = new types::Double(pPolyY0->getCoef()->getCols(), pPolyY0->getCoef()->getRows());
205 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ode", iPos + 1);
206 return types::Function::Error;
211 if (in[iPos]->isDouble() == false)
213 Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "ode", iPos + 1);
214 return types::Function::Error;
217 pDblT0 = in[iPos]->getAs<types::Double>();
219 if (pDblT0->isScalar() == false)
221 Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "ode", iPos + 1);
222 return types::Function::Error;
227 if (in[iPos]->isDouble() == false)
229 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ode", iPos + 1);
230 return types::Function::Error;
233 pDblT = in[iPos]->getAs<types::Double>();
236 DifferentialEquationFunctions* deFunctionsManager = new DifferentialEquationFunctions(L"ode");
237 DifferentialEquation::addDifferentialEquationFunctions(deFunctionsManager);
238 deFunctionsManager->setOdeYRows(pDblY0->getRows());
239 deFunctionsManager->setOdeYCols(pDblY0->getCols());
241 YSize = (int*)MALLOC(sizeOfYSize * sizeof(int));
242 *YSize = pDblY0->getSize();
243 pdYData = (double*)MALLOC(pDblY0->getSize() * sizeof(double));
244 C2F(dcopy)(YSize, pDblY0->get(), &one, pdYData, &one);
250 Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ode", 5);
251 DifferentialEquation::removeDifferentialEquationFunctions();
254 return types::Function::Error;
257 if (in[4]->isCallable() == false && in[4]->isString() == false && in[4]->isList() == false)
259 Scierror(999, _("%s: Wrong type for input argument #%d: A function expected.\n"), "ode", 5);
260 DifferentialEquation::removeDifferentialEquationFunctions();
263 return types::Function::Error;
267 for (iPos++; iPos < (int)in.size(); iPos++)
269 if (in[iPos]->isDouble())
271 if (pDblRtol == NULL && bFuncF == false)
273 pDblRtol = in[iPos]->getAs<types::Double>();
274 if (pDblRtol->getSize() != pDblY0->getSize() && pDblRtol->isScalar() == false)
276 Scierror(267, _("%s: Arg %d and arg %d must have equal dimensions.\n"), "ode", pStrType ? 2 : 1, iPos + 1);
277 DifferentialEquation::removeDifferentialEquationFunctions();
280 return types::Function::Error;
283 else if (pDblAtol == NULL && bFuncF == false)
285 pDblAtol = in[iPos]->getAs<types::Double>();
286 if (pDblAtol->getSize() != pDblY0->getSize() && pDblAtol->isScalar() == false)
288 Scierror(267, _("%s: Arg %d and arg %d must have equal dimensions.\n"), "ode", pStrType ? 2 : 1, iPos + 1);
289 DifferentialEquation::removeDifferentialEquationFunctions();
292 return types::Function::Error;
295 else if (pDblNg == NULL && bFuncF == true && meth == 3)
297 pDblNg = in[iPos]->getAs<types::Double>();
299 else if (pDblW == NULL && bFuncF == true && (bFuncG == true || meth != 3))
301 if ((int)in.size() == iPos + 2)
303 if (in[iPos + 1]->isDouble() == false)
305 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ode", iPos + 2);
306 DifferentialEquation::removeDifferentialEquationFunctions();
309 return types::Function::Error;
312 pDblW = in[iPos]->getAs<types::Double>();
313 pDblIw = in[iPos + 1]->getAs<types::Double>();
318 Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ode", iPos + 2);
319 DifferentialEquation::removeDifferentialEquationFunctions();
322 return types::Function::Error;
327 Scierror(999, _("%s: Wrong type for input argument #%d: A function expected.\n"), "ode", iPos + 1);
328 DifferentialEquation::removeDifferentialEquationFunctions();
331 return types::Function::Error;
334 else if (in[iPos]->isCallable())
336 types::Callable* pCall = in[iPos]->getAs<types::Callable>();
339 deFunctionsManager->setFFunction(pCall);
342 else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
344 deFunctionsManager->setJacFunction(pCall);
347 else if (bFuncG == false && meth == 3)
349 deFunctionsManager->setGFunction(pCall);
354 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ode", iPos + 1);
355 DifferentialEquation::removeDifferentialEquationFunctions();
358 return types::Function::Error;
361 else if (in[iPos]->isString())
363 types::String* pStr = in[iPos]->getAs<types::String>();
368 bOK = deFunctionsManager->setFFunction(pStr);
371 else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
373 bOK = deFunctionsManager->setJacFunction(pStr);
376 else if (bFuncG == false && meth == 3)
378 bOK = deFunctionsManager->setGFunction(pStr);
383 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ode", iPos + 1);
384 DifferentialEquation::removeDifferentialEquationFunctions();
387 return types::Function::Error;
392 char* pst = wide_string_to_UTF8(pStr->get(0));
393 Scierror(50, _("%s: Subroutine not found: %s\n"), "ode", pst);
395 DifferentialEquation::removeDifferentialEquationFunctions();
398 return types::Function::Error;
401 else if (in[iPos]->isList())
403 types::List* pList = in[iPos]->getAs<types::List>();
405 if (pList->getSize() == 0)
407 Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "ode", iPos + 1, "(string empty)");
408 DifferentialEquation::removeDifferentialEquationFunctions();
411 return types::Function::Error;
414 if (bFuncF && (bFuncJac || pDblNg) && (bFuncG || meth != 3))
416 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ode", iPos + 1);
417 DifferentialEquation::removeDifferentialEquationFunctions();
420 return types::Function::Error;
423 if (pList->get(0)->isString())
425 types::String* pStr = pList->get(0)->getAs<types::String>();
431 bOK = deFunctionsManager->setFFunction(pStr);
432 sizeOfpdYData = *YSize;
434 else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
437 bOK = deFunctionsManager->setJacFunction(pStr);
438 if (sizeOfpdYData == 0)
440 sizeOfpdYData = *YSize;
443 else if (bFuncG == false && meth == 3)
446 bOK = deFunctionsManager->setGFunction(pStr);
447 if (sizeOfpdYData == 0)
449 sizeOfpdYData = *YSize;
455 char* pst = wide_string_to_UTF8(pStr->get(0));
456 Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "ode", iPos + 1, pst);
458 DifferentialEquation::removeDifferentialEquationFunctions();
461 return types::Function::Error;
464 int* sizeTemp = YSize;
465 int totalSize = sizeOfpdYData;
467 YSize = (int*)MALLOC((sizeOfYSize + pList->getSize() - 1) * sizeof(int));
468 memcpy(YSize, sizeTemp, sizeOfYSize * sizeof(int));
470 std::vector<types::Double*> vpDbl;
471 for (int iter = 0; iter < pList->getSize() - 1; iter++)
473 if (pList->get(iter + 1)->isDouble() == false)
475 Scierror(999, _("%s: Wrong type for input argument #%d: Argument %d in the list must be a matrix.\n"), "ode", iPos + 1, iter + 1);
476 DifferentialEquation::removeDifferentialEquationFunctions();
479 return types::Function::Error;
482 vpDbl.push_back(pList->get(iter + 1)->getAs<types::Double>());
483 YSize[sizeOfYSize + iter] = vpDbl[iter]->getSize();
484 totalSize += YSize[sizeOfYSize + iter];
487 double* pdYDataTemp = pdYData;
488 pdYData = (double*)MALLOC(totalSize * sizeof(double));
489 C2F(dcopy)(&sizeOfpdYData, pdYDataTemp, &one, pdYData, &one);
491 int position = sizeOfpdYData;
492 for (int iter = 0; iter < pList->getSize() - 1; iter++)
494 C2F(dcopy)(&YSize[sizeOfYSize + iter], vpDbl[iter]->get(), &one, &pdYData[position], &one);
495 position += vpDbl[iter]->getSize();
498 sizeOfpdYData = totalSize;
499 sizeOfYSize += pList->getSize() - 1;
503 else if (pList->get(0)->isCallable())
508 deFunctionsManager->setFFunction(pList->get(0)->getAs<types::Callable>());
509 for (int iter = 1; iter < pList->getSize(); iter++)
511 deFunctionsManager->setFArgs(pList->get(iter)->getAs<types::InternalType>());
514 else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
517 deFunctionsManager->setJacFunction(pList->get(0)->getAs<types::Callable>());
518 for (int iter = 1; iter < pList->getSize(); iter++)
520 deFunctionsManager->setJacArgs(pList->get(iter)->getAs<types::InternalType>());
523 else if (bFuncG == false && meth == 3)
526 deFunctionsManager->setGFunction(pList->get(0)->getAs<types::Callable>());
527 for (int iter = 1; iter < pList->getSize(); iter++)
529 deFunctionsManager->setGArgs(pList->get(iter)->getAs<types::InternalType>());
535 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);
536 DifferentialEquation::removeDifferentialEquationFunctions();
539 return types::Function::Error;
544 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a function expected.\n"), "ode", iPos + 1);
545 DifferentialEquation::removeDifferentialEquationFunctions();
548 return types::Function::Error;
554 int val = (meth == 3) ? 3 : 1;
555 Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ode", in.size() + val);
556 DifferentialEquation::removeDifferentialEquationFunctions();
559 return types::Function::Error;
561 if (pDblNg == NULL && meth == 3)
563 Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ode", in.size() + 2);
564 DifferentialEquation::removeDifferentialEquationFunctions();
567 return types::Function::Error;
569 if (bFuncG == false && meth == 3)
571 Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ode", in.size() + 1);
572 DifferentialEquation::removeDifferentialEquationFunctions();
575 return types::Function::Error;
578 // *** Initialization. ***
581 int itask = 1; // itask = 1 for normal computation of output values of y at t = tout.
582 int istate = 1; // istate = integer flag (input and output). set istate = 1.
583 int iopt = 0; // iopt = 0 to indicate no optional inputs used.
584 int jt = 2; // jacobian type indicator. set jt = 2.
589 double* rwork = NULL;
594 // contain ls0001, lsa001 and eh0001 structures
595 double* dStructTab = NULL;
596 int* iStructTab = NULL;
597 int dStructTabSize = 0;
598 int iStructTabSize = 0;
600 int rwSize = 0; // rwSize = dStructTab + rworkSize
601 int iwSize = 0; // iwSize = iStructTab + iworkSize
603 // structures used by lsoda and lsode
604 double* ls0001d = &(C2F(ls0001).tret);
605 int* ls0001i = &(C2F(ls0001).illin);
606 double* lsa001d = &(C2F(lsa001).tsw);
607 int* lsa001i = &(C2F(lsa001).insufr);
608 double* lsr001d = &(C2F(lsr001).rownr3[0]);
609 int* lsr001i = &(C2F(lsr001).lg0);
610 int* eh0001i = &(C2F(eh0001).mesflg);
613 types::InternalType* pIT = symbol::Context::getInstance()->get(symbol::Symbol(L"%ODEOPTIONS"));
614 if (pIT != NULL && pIT->isDouble())
616 pDblOdeOptions = pIT->getAs<types::Double>();
617 if (pDblOdeOptions->getSize() == 12)
620 itask = (int)pDblOdeOptions->get(0);
621 jt = (int)pDblOdeOptions->get(5);
622 ml = (int)pDblOdeOptions->get(10);
623 mu = (int)pDblOdeOptions->get(11);
627 if (iopt == 1 && (pDblOdeOptions->get(4) > pDblOdeOptions->get(3))) // hmin > hmax ?
629 Scierror(9999, _("%s: Wrong value of hmin and hmax: hmin = %d is greater than hmax = %d.\n"), "ode", pDblOdeOptions->get(4), pDblOdeOptions->get(3));
630 DifferentialEquation::removeDifferentialEquationFunctions();
633 return types::Function::Error;
636 if (jt < 0 || jt > 5)
638 Scierror(9999, _("%s: Wrong value of Jacobian type: A number between %d and %d expected.\n"), "ode", 0, 5);
639 DifferentialEquation::removeDifferentialEquationFunctions();
642 return types::Function::Error;
645 if (iopt == 0 && bFuncJac)
650 if (bFuncJac && (jt == 2 || jt == 5) && getWarningMode())
652 sciprint(_("%s: Warning: Jacobian is given, but not used.\n"), "ode");
655 if (bFuncJac == false && (jt == 1 || jt == 4))
657 if (getWarningMode())
659 sciprint(_("%s: Warning: No Jacobian external given, but one is required by %ODEOPTIONS(6) value. Jacobian will be estimated.\n"), "ode");
665 //compute itol and set the tolerances rtol and atol.
671 if (pDblRtol->isScalar())
673 rtol = (double*)MALLOC(sizeof(double));
674 *rtol = pDblRtol->get(0);
678 rtol = pDblRtol->get();
684 rtol = (double*)MALLOC(sizeof(double));
685 if (meth == 6 || meth == 7)
697 if (pDblAtol->isScalar())
699 atol = (double*)MALLOC(sizeof(double));
700 *atol = pDblAtol->get(0);
704 atol = pDblAtol->get();
710 atol = (double*)MALLOC(sizeof(double));
711 if (meth == 6 || meth == 7)
721 // Compute rwork, iwork size.
724 if (pDblW) // structure ls0001 have been restored.
726 nyh = C2F(ls0001).nyh;
734 mxordn = (int)pDblOdeOptions->get(7);
735 mxords = (int)pDblOdeOptions->get(8);
738 if (mxordn > 12 || mxords > 5 || mxordn < 1 || mxords < 1)
740 if (getWarningMode())
742 sciprint(_("%s: Warning: Wrong value for maximun 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);
756 case 3: // lsodar (root)
758 // lrn = 20 + nyh*(mxordn+1) + 3*neq + 3*ng,
759 // lrs = 20 + nyh*(mxords+1) + 3*neq + lmat + 3*ng,
761 // nyh = the initial value of neq,
762 // mxordn = 12, unless a smaller value is given as an
764 // mxords = 5, unless a smaller value is given as an
766 // lmat = length of matrix work space..
767 // lmat = neq**2 + 2 if jt = 1 or 2,
768 // lmat = (2*ml + mu + 1)*neq + 2 if jt = 4 or 5.
770 lrn = 3 * (int)pDblNg->get(0);
772 dStructTabSize = 246 - 241;
773 iStructTabSize = 59 - 50;
774 jroot = (int*)MALLOC((int)pDblNg->get(0) * sizeof(int));
778 // lrn = 20 + nyh*(mxordn+1) + 3*neq,
779 // lrs = 20 + nyh*(mxords+1) + 3*neq + lmat,
781 // nyh = the initial value of neq,
782 // mxordn = 12, unless a smaller value is given as an
784 // mxords = 5, unless a smaller value is given as an
786 // lmat = length of matrix work space..
787 // lmat = neq**2 + 2 if jt = 1 or 2,
788 // lmat = (2*ml + mu + 1)*neq + 2 if jt = 4 or 5.
790 if (jt == 1 || jt == 2)
792 lmat = (*YSize) * (*YSize) + 2; // if jt = 1 or 2,
794 else if (jt == 4 || jt == 5)
796 //ml and mu = -1 in all cases
797 lmat = (2 * ml + mu + 1) * (*YSize) + 2; // if jt = 4 or 5.
800 lrn += 20 + nyh * (mxordn + 1) + 3 * (*YSize);
801 lrs += 20 + nyh * (mxords + 1) + 3 * (*YSize) + lmat;
803 rworkSize = max(lrn, lrs);
804 iworkSize = 20 + *YSize;
806 dStructTabSize += 241;
807 iStructTabSize += 50;
811 case 1: // lsode (adams) (non stiff)
815 case 2: // lsode (stiff)
817 // 20 + nyh*(maxord + 1) + 3*neq + lmat
819 // nyh = the initial value of neq,
820 // maxord = 12 (if meth = 1) or 5 (if meth = 2) (unless a
821 // smaller value is given as an optional input),
822 // lmat = 0 if miter = 0,
823 // lmat = neq**2 + 2 if miter is 1 or 2,
824 // lmat = neq + 2 if miter = 3, and
825 // lmat = (2*ml+mu+1)*neq + 2 if miter is 4 or 5.
827 if (jt == 1 || jt == 2)
829 lmat = (*YSize) * (*YSize) + 2;
835 else if (jt == 4 || jt == 5)
837 lmat = (2 * ml + mu + 1) * (*YSize) + 2;
840 rworkSize = 20 + nyh * (maxord + 1) + 3 * (*YSize) + lmat;
841 iworkSize = 20; // if jt = 0 or 3
842 if (jt == 1 || jt == 2 || jt == 4 || jt == 5) // iSize = 20 + neq otherwise
844 iworkSize += (*YSize);
847 dStructTabSize = 219;
851 case 4: // lsdisc (discrete)
857 case 5: // lsrgk (rk)
859 rworkSize = 9 * (*YSize);
863 case 6: // rkf45 (rkf)
865 rworkSize = 3 + 8 * (*YSize);
869 case 7: // rksimp (fix)
871 rworkSize = 3 + 8 * (*YSize);
877 rwSize = rworkSize + dStructTabSize;
878 iwSize = iworkSize + iStructTabSize;
880 rwork = (double*)MALLOC(rworkSize * sizeof(double));
881 iwork = (int*)MALLOC(iworkSize * sizeof(int));
887 if (pDblW->getSize() != rwSize || pDblIw->getSize() != iwSize)
889 Scierror(9999, _("%s: Wrong size for w and iw: w = %d and iw = %d expected.\n"), "ode", rwSize, iwSize);
890 DifferentialEquation::removeDifferentialEquationFunctions();
899 if (itol == 1 || itol == 3)
907 return types::Function::Error;
910 istate = 2; // 1 means this is the first call | 2 means this is not the first call
912 dStructTab = (double*)MALLOC(dStructTabSize * sizeof(double));
913 iStructTab = (int*)MALLOC(iStructTabSize * sizeof(int));
915 C2F(dcopy)(&rworkSize, pDblW->get(), &one, rwork, &one);
916 C2F(dcopy)(&dStructTabSize, pDblW->get() + rworkSize, &one, dStructTab, &one);
918 for (int i = 0; i < iworkSize; i++)
920 iwork[i] = (int)pDblIw->get(i);
923 for (int i = 0; i < iStructTabSize; i++)
925 iStructTab[i] = (int)pDblIw->get(i + iworkSize);
930 //if iopt = 0 default value are used without set it.
935 rwork[0] = pDblOdeOptions->get(1); // tcrit
938 rwork[4] = pDblOdeOptions->get(2); // h0 =0
939 rwork[5] = pDblOdeOptions->get(3); // hmax = %inf
940 rwork[6] = pDblOdeOptions->get(4); // hmin = 0
942 if (jt == 4 || jt == 5)
944 iwork[0] = (int)pDblOdeOptions->get(10);// ml = -1
945 iwork[1] = (int)pDblOdeOptions->get(11);// mu = -1
948 if (meth == 0 || meth == 3)
951 iwork[4] = (int)pDblOdeOptions->get(9); // ixpr = 0
952 iwork[7] = (int)pDblOdeOptions->get(7); // mxordn = 12
953 iwork[8] = (int)pDblOdeOptions->get(8); // mxords = 5
960 iwork[4] = (int)pDblOdeOptions->get(7); // mxordn = 12
964 iwork[4] = (int)pDblOdeOptions->get(8); // mxords = 5
968 iwork[5] = (int)pDblOdeOptions->get(6); // mxstep = 500
969 iwork[6] = 10; // mxhnil = 10 maximum number of messages printed per problem
974 if (pDblOdeOptions && pDblOdeOptions->get(9) == 1)
976 sciprint(_("itask = %d\tmeth = %d\tjactyp = %d\tml = %d\tmu = %d\tiopt = %d\n"), itask, meth, jt, ml, mu, iopt);
977 sciprint(_("tcrit = %lf\th0 = %lf\thmax = %lf\thmin = %lf\n"), pDblOdeOptions->get(1), pDblOdeOptions->get(2), pDblOdeOptions->get(3), pDblOdeOptions->get(4));
988 //restore ls0001 and eh0001 from w (rwork) and iw (iwork).
989 C2F(dcopy)(&dSize, dStructTab, &one, ls0001d, &one);
990 memcpy(ls0001i, iStructTab, 39 * sizeof(int));
995 //restore lsa001 from w (rwork) and iw (iwork).
996 if (meth == 0 || meth == 3)
999 C2F(dcopy)(&dSize, dStructTab + dPos, &one, lsa001d, &one);
1000 memcpy(lsa001i, &iStructTab[iPos], 9 * sizeof(int));
1010 C2F(dcopy)(&dSize, dStructTab + dPos, &one, lsr001d, &one);
1011 memcpy(lsr001i, &iStructTab[iPos], 9 * sizeof(int));
1016 // *** Perform operation. ***
1018 double t0 = pDblT0->get(0);
1019 bool bOneStep = false;
1021 if (itask == 5 || itask == 3 || itask == 2)
1024 if (getWarningMode() && pDblT->isScalar() == false)
1026 sciprint(_("itask = %d: At most one value of t is allowed, the last element of t is used.\n"), itask);
1032 std::list<double*> pDblYOutList = std::list<double*>();
1033 std::list<double> pDblTOutList = std::list<double>();
1034 int iLastT = pDblT->getSize() - 1;
1035 double t = pDblT->get(iLastT);
1036 int iDir = t - t0 < 0 ? -1 : 1;
1040 std::string strMeth;
1045 case 1: // lsode (adams)
1046 case 2: // lsode (stiff)
1049 int jacType = 10 * meth + jt;
1050 C2F(lsode)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jacType);
1056 int ng = (int)pDblNg->get(0);
1057 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);
1060 case 4: // lsdisc (discrete)
1063 C2F(lsdisc)(ode_f, YSize, pdYData, &t0, &t, rwork, &rworkSize, &istate);
1066 case 5: // lsrgk (rk)
1069 C2F(lsrgk)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1072 case 6: // rkf45 (rkf)
1075 C2F(rkf45)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1078 case 7: // rksimp (fix)
1081 C2F(rksimp)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1084 default: // case 0: // lsoda
1087 C2F(lsoda)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt);
1092 err = checkOdeError(meth, istate);
1095 Scierror(999, _("%s: %s exit with state %d.\n"), "ode", strMeth.c_str(), istate);
1098 catch (ast::ScilabError &e)
1100 char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
1101 sciprint(_("%s: exception caught in '%s' subroutine.\n"), "ode", strMeth.c_str());
1102 Scierror(999, pstrMsg);
1107 // FREE allocated data
1108 if (err == 1) // error case
1110 DifferentialEquation::removeDifferentialEquationFunctions();
1127 if (itol == 1 || itol == 3)
1135 return types::Function::Error;
1138 pDblYOutList.push_back(pdYData);
1139 pDblTOutList.push_back(t0);
1141 if (err == 2) // warning case
1143 if (getWarningMode())
1145 sciprint(_("Integration was stoped at t = %lf.\n"), t0);
1150 if (meth == 3 && istate == 3)
1152 // istate == 3 means the integration was successful, and one or more
1153 // roots were found before satisfying the stop condition
1154 // specified by itask. see jroot.
1155 if (getWarningMode())
1157 sciprint(_("%s: Warning: At t = %lf, y is a root, jroot = "), "ode", t0);
1158 for (int k = 0; k < pDblNg->get(0); k++)
1160 sciprint("\t%d", jroot[k]);
1168 while ((t0 - t) * iDir < 0);
1170 int iSizeList = (int)pDblYOutList.size();
1172 if(pPolyY0) // pPolyY0 is scalar
1174 // result format = [2 x iSizeList]
1175 // first lime is the time of result stored in second line.
1176 int size = 2 * iSizeList;
1177 int* iRanks = (int*)MALLOC(size * sizeof(int));
1178 int iMaxRank = pPolyY0->getMaxRank();
1179 for(int i = 0; i < iSizeList; i++)
1181 iRanks[i * 2] = 1; // time rank
1182 iRanks[i * 2 + 1] = iMaxRank; // result rank
1185 pPolyYOut = new types::Polynom(pPolyY0->getVariableName(), 2, iSizeList, iRanks);
1186 pDblYOut = new types::Double(pDblY0->getRows(), 1);
1187 pDblTOut = new types::Double(1, 1);
1189 for(int i = 0; i < iSizeList; i++)
1191 // pDblY0 contain pPolyY0 (ie: pPolyY0 = 2 + x => pDblY0 = [2 1])
1192 for(int j = 0; j < pDblY0->getRows(); j++)
1194 pDblYOut->set(j, pDblYOutList.front()[j]);
1196 pDblTOut->set(0, pDblTOutList.front());
1198 pPolyYOut->setCoef(i * 2, pDblTOut);
1199 pPolyYOut->setCoef(i * 2 + 1, pDblYOut);
1201 pDblYOutList.pop_front();
1202 pDblTOutList.pop_front();
1208 pDblYOut = new types::Double(pDblY0->getRows() + 1, iSizeList);
1209 for (int i = 0; i < iSizeList; i++)
1211 pDblYOut->set(i * (*YSize + 1), pDblTOutList.front());
1212 for (int j = 0; j < *YSize; j++)
1214 pDblYOut->set(i * (*YSize + 1) + (j + 1), pDblYOutList.front()[j]);
1216 pDblYOutList.pop_front();
1217 pDblTOutList.pop_front();
1226 int size = pDblT->getSize();
1227 int* iRanks = (int*)MALLOC(size * sizeof(int));
1228 int iMaxRank = pPolyY0->getMaxRank();
1229 for(int i = 0; i < size; i++)
1231 iRanks[i] = iMaxRank;
1234 pPolyYOut = new types::Polynom(pPolyY0->getVariableName(), 1, pDblT->getSize(), iRanks);
1236 pDblYOut = new types::Double(pDblY0->getRows(), 1);
1241 pDblYOut = new types::Double(pDblY0->getRows(), pDblY0->getCols() * pDblT->getSize());
1243 bool bBreak = false;
1244 for (int i = 0; i < pDblT->getSize(); i++)
1246 double t = pDblT->get(i);
1247 std::string strMeth;
1249 if (itask >= 4 && t > rwork[0]) // rwork[0] => tcrit
1259 case 1: // lsode (adams)
1260 case 2: // lsode (stiff)
1263 int jacType = 10 * meth + jt;
1264 C2F(lsode)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jacType);
1270 int ng = (int)pDblNg->get(0);
1271 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);
1274 case 4: // lsdisc (discrete)
1277 C2F(lsdisc)(ode_f, YSize, pdYData, &t0, &t, rwork, &rworkSize, &istate);
1280 case 5: // lsrgk (rk)
1283 C2F(lsrgk)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1286 case 6: // rkf45 (rkf)
1289 C2F(rkf45)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1292 case 7: // rksimp (fix)
1295 C2F(rksimp)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1298 default: // case 0: // lsoda
1301 C2F(lsoda)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt);
1306 err = checkOdeError(meth, istate);
1309 Scierror(999, _("%s: %s exit with state %d.\n"), "ode", strMeth.c_str(), istate);
1312 catch (ast::ScilabError &e)
1314 char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
1315 sciprint(_("%s: exception caught in '%s' subroutine.\n"), "ode", strMeth.c_str());
1316 Scierror(999, pstrMsg);
1321 // FREE allocated data
1322 if (err == 1) // error case
1324 DifferentialEquation::removeDifferentialEquationFunctions();
1341 if (itol == 1 || itol == 3)
1349 return types::Function::Error;
1354 for(int j = 0; j < pDblY0->getRows(); j++)
1356 pDblYOut->set(j, pdYData[j]);
1358 pPolyYOut->setCoef(i, pDblYOut);
1363 for (int j = 0; j < *YSize; j++)
1365 pDblYOut->set(i * (*YSize) + j, pdYData[j]);
1369 if (err == 2) // warning case
1371 if (getWarningMode())
1373 sciprint(_("Integration was stoped at t = %lf.\n"), t0);
1376 types::Double* pDblYOutTemp = pDblYOut;
1377 pDblYOut = new types::Double(pDblYOutTemp->getRows(), i + 1);
1378 for (int k = 0; k < pDblYOut->getSize(); k++)
1380 pDblYOut->set(k, pDblYOutTemp->get(k));
1382 delete pDblYOutTemp;
1387 if (meth == 3 && istate == 3)
1389 // istate == 3 means the integration was successful, and one or more
1390 // roots were found before satisfying the stop condition
1391 // specified by itask. see jroot.
1393 if (getWarningMode())
1395 sciprint(_("%s: Warning: At t = %lf, y is a root, jroot = "), "ode", t0);
1396 for (int k = 0; k < pDblNg->get(0); k++)
1398 sciprint("\t%d", jroot[k]);
1403 types::Double* pDblYOutTemp = pDblYOut;
1404 pDblYOut = new types::Double(pDblYOutTemp->getRows(), i + 1);
1405 for (int k = 0; k < pDblYOut->getSize(); k++)
1407 pDblYOut->set(k, pDblYOutTemp->get(k));
1409 delete pDblYOutTemp;
1413 if (itask >= 4 && bBreak)
1415 types::Double* pDblYOutTemp = pDblYOut;
1416 pDblYOut = new types::Double(pDblYOutTemp->getRows(), i + 1);
1417 for (int k = 0; k < pDblYOut->getSize(); k++)
1419 pDblYOut->set(k, pDblYOutTemp->get(k));
1421 delete pDblYOutTemp;
1429 if (_iRetCount > 2) //save ls0001 and eh0001 following w and iw.
1435 if (dStructTab == NULL)
1437 dStructTab = (double*)MALLOC(dStructTabSize * sizeof(double));
1438 iStructTab = (int*)MALLOC(iStructTabSize * sizeof(int));
1441 C2F(dcopy)(&dSize, ls0001d, &one, dStructTab, &one);
1442 memcpy(iStructTab, ls0001i, 39 * sizeof(int));
1448 if (meth == 0 || meth == 3)
1451 C2F(dcopy)(&dSize, lsa001d, &one, dStructTab + dPos, &one);
1452 memcpy(&iStructTab[iPos], lsa001i, 9 * sizeof(int));
1462 C2F(dcopy)(&dSize, lsr001d, &one, dStructTab + dPos, &one);
1463 memcpy(&iStructTab[iPos], lsr001i, 9 * sizeof(int));
1468 memcpy(&iStructTab[iPos], eh0001i, 2 * sizeof(int));
1472 // *** Return result in Scilab. ***
1476 out.push_back(pPolyYOut); // y
1481 out.push_back(pDblYOut); // y
1484 if (meth == 3 && _iRetCount >= 2)
1489 for (int i = 0; i < pDblNg->get(0); i++)
1497 types::Double* pDblRd = new types::Double(1, sizeOfRd);
1498 //rd: The first entry contains the stopping time.
1499 pDblRd->set(0, C2F(lsr001).tlast);
1500 for (int i = 0; i < pDblNg->get(0); i++)
1505 pDblRd->set(k, (double)i + 1);
1508 out.push_back(pDblRd); // rd
1511 if ((meth < 3 && _iRetCount == 3) || (meth == 3 && _iRetCount == 4))
1513 types::Double* pDblWOut = new types::Double(1, rwSize);
1514 C2F(dcopy)(&rworkSize, rwork, &one, pDblWOut->get(), &one);
1515 C2F(dcopy)(&dStructTabSize, dStructTab, &one, pDblWOut->get() + rworkSize, &one);
1517 types::Double* pDblIwOut = new types::Double(1, iwSize);
1518 for (int i = 0; i < iworkSize; i++)
1520 pDblIwOut->set(i, (double)iwork[i]);
1523 for (int i = 0; i < iStructTabSize; i++)
1525 pDblIwOut->set(iworkSize + i, (double)iStructTab[i]);
1528 out.push_back(pDblWOut); // w
1529 out.push_back(pDblIwOut); // iw
1533 if (itol == 1 || itol == 3) // atol is scalar
1538 if (itol < 3) // rtol is scalar
1559 DifferentialEquation::removeDifferentialEquationFunctions();
1561 return types::Function::OK;