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 "sci_malloc.h"
28 #include "localization.h"
30 #include "elem_common.h"
31 #include "configvariable_interface.h"
33 #include "common_structure.h"
34 #include "scifunctions.h"
37 /*--------------------------------------------------------------------------*/
38 types::Function::ReturnValue sci_odedc(types::typed_list &in, int _iRetCount, types::typed_list &out)
41 types::String* pStrType = NULL;
42 const wchar_t* wcsType = L"lsoda";
46 types::Double* pDblY0 = NULL;
47 double* pdYData = NULL; // contain y0 following by all args data in list case.
48 int sizeOfpdYData = 0;
51 types::Double* pDblT0 = NULL;
52 types::Double* pDblT = NULL;
53 types::Double* pDblRtol = NULL;
54 types::Double* pDblAtol = NULL;
55 types::Double* pDblNg = NULL;
56 types::Double* pDblW = NULL;
57 types::Double* pDblIw = NULL;
58 types::Double* pDblNd = NULL;
59 types::Double* pDblStdel = NULL;
62 types::Double* pDblOdeOptions = NULL;
65 types::Double* pDblYOut = NULL;
67 // Indicate if the function is given.
69 bool bFuncJac = false;
72 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"), "odedc", 6);
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"), "odedc", 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"), "odedc", 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"), "odedc", 1, 2, 4);
154 return types::Function::Error;
161 Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "odedc", 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"), "odedc", iPos + 1);
174 return types::Function::Error;
177 if (pDblY0->getCols() != 1)
179 Scierror(999, _("%s: Wrong size for input argument #%d: A real colunm vector expected (n x 1).\n"), "odedc", iPos + 1);
180 return types::Function::Error;
185 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix\n"), "odedc", iPos + 1);
186 return types::Function::Error;
191 if (in[iPos]->isDouble() == false)
193 Scierror(999, _("%s: Wrong type for input argument #%d: A real scalar expected.\n"), "odedc", iPos + 1);
194 return types::Function::Error;
197 pDblNd = in[iPos]->getAs<types::Double>();
199 if (pDblNd->isComplex() || pDblNd->isScalar() == false)
201 Scierror(999, _("%s: Wrong type for input argument #%d: A real scalar expected.\n"), "odedc", iPos + 1);
202 return types::Function::Error;
205 if (pDblNd->get(0) > pDblY0->getSize())
207 Scierror(999, _("%s: Wrong value for input argument #%d: Value must not exceeds dimension of argument %d.\n"), "odedc", iPos + 1, iPos);
208 return types::Function::Error;
211 // stdel => hstep and/or delta with delta=0 as default value
213 if (in[iPos]->isDouble() == false)
215 Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "odedc", iPos + 1);
216 return types::Function::Error;
219 pDblStdel = in[iPos]->getAs<types::Double>();
221 if (pDblStdel->isComplex())
223 Scierror(999, _("%s: Wrong type for input argument #%d: A real scalar expected.\n"), "odedc", iPos + 1);
224 return types::Function::Error;
227 if (pDblStdel->getSize() > 2)
229 Scierror(999, _("%s: Wrong size for input argument #%d: %d or %d values expected.\n"), "odedc", iPos + 1, 1, 2);
230 return types::Function::Error;
235 if (in[iPos]->isDouble() == false)
237 Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "odedc", iPos + 1);
238 return types::Function::Error;
241 pDblT0 = in[iPos]->getAs<types::Double>();
243 if (pDblT0->isScalar() == false)
245 Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "odedc", iPos + 1);
246 return types::Function::Error;
251 if (in[iPos]->isDouble() == false)
253 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "odedc", iPos + 1);
254 return types::Function::Error;
257 pDblT = in[iPos]->getAs<types::Double>();
259 // size of continuous part
260 sizeYc = pDblY0->getSize() - (int)pDblNd->get(0);
263 Scierror(999, _("%s: Wrong value for input argument #%d: Value of nd exceeds dimension of y0.\n"), "odedc", 2);
264 return types::Function::Error;
268 DifferentialEquationFunctions* deFunctionsManager = new DifferentialEquationFunctions(L"odedc");
269 DifferentialEquation::addDifferentialEquationFunctions(deFunctionsManager);
270 deFunctionsManager->setOdedcYDSize((int)pDblNd->get(0));
271 deFunctionsManager->setOdeYRows(pDblY0->getRows());
272 deFunctionsManager->setOdeYCols(pDblY0->getCols());
274 YSize = (int*)MALLOC(sizeOfYSize * sizeof(int));
275 *YSize = pDblY0->getSize();
276 pdYData = (double*)MALLOC(pDblY0->getSize() * sizeof(double));
277 C2F(dcopy)(YSize, pDblY0->get(), &one, pdYData, &one);
283 Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "odedc", 7);
284 DifferentialEquation::removeDifferentialEquationFunctions();
287 return types::Function::Error;
290 if (in[6]->isCallable() == false && in[6]->isString() == false && in[6]->isList() == false)
292 Scierror(999, _("%s: Wrong type for input argument #%d: A function expected.\n"), "odedc", 7);
293 DifferentialEquation::removeDifferentialEquationFunctions();
296 return types::Function::Error;
300 for (iPos++; iPos < in.size(); iPos++)
302 if (in[iPos]->isDouble())
304 if (pDblRtol == NULL && bFuncF == false)
306 pDblRtol = in[iPos]->getAs<types::Double>();
307 if (pDblRtol->getSize() != pDblY0->getSize() && pDblRtol->isScalar() == false)
309 Scierror(267, _("%s: Arg %d and arg %d must have equal dimensions.\n"), "odedc", pStrType ? 2 : 1, iPos + 1);
310 DifferentialEquation::removeDifferentialEquationFunctions();
313 return types::Function::Error;
316 else if (pDblAtol == NULL && bFuncF == false)
318 pDblAtol = in[iPos]->getAs<types::Double>();
319 if (pDblAtol->getSize() != pDblY0->getSize() && pDblAtol->isScalar() == false)
321 Scierror(267, _("%s: Arg %d and arg %d must have equal dimensions.\n"), "odedc", pStrType ? 2 : 1, iPos + 1);
322 DifferentialEquation::removeDifferentialEquationFunctions();
325 return types::Function::Error;
328 else if (pDblNg == NULL && bFuncF == true && meth == 3)
330 pDblNg = in[iPos]->getAs<types::Double>();
332 else if (pDblW == NULL && bFuncF == true && (bFuncG == true || meth != 3))
334 if (in.size() == iPos + 2)
336 if (in[iPos + 1]->isDouble() == false)
338 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "odedc", iPos + 2);
339 DifferentialEquation::removeDifferentialEquationFunctions();
342 return types::Function::Error;
345 pDblW = in[iPos]->getAs<types::Double>();
346 pDblIw = in[iPos + 1]->getAs<types::Double>();
351 Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "odedc", iPos + 2);
352 DifferentialEquation::removeDifferentialEquationFunctions();
355 return types::Function::Error;
360 Scierror(999, _("%s: Wrong type for input argument #%d: A function expected.\n"), "odedc", iPos + 1);
361 DifferentialEquation::removeDifferentialEquationFunctions();
364 return types::Function::Error;
367 else if (in[iPos]->isCallable())
369 types::Callable* pCall = in[iPos]->getAs<types::Callable>();
372 deFunctionsManager->setFFunction(pCall);
375 else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
377 deFunctionsManager->setJacFunction(pCall);
380 else if (bFuncG == false && meth == 3)
382 deFunctionsManager->setGFunction(pCall);
387 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "odedc", iPos + 1);
388 DifferentialEquation::removeDifferentialEquationFunctions();
391 return types::Function::Error;
394 else if (in[iPos]->isString())
396 types::String* pStr = in[iPos]->getAs<types::String>();
401 bOK = deFunctionsManager->setFFunction(pStr);
404 else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
406 bOK = deFunctionsManager->setJacFunction(pStr);
409 else if (bFuncG == false && meth == 3)
411 bOK = deFunctionsManager->setGFunction(pStr);
416 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "odedc", iPos + 1);
417 DifferentialEquation::removeDifferentialEquationFunctions();
420 return types::Function::Error;
425 char* pst = wide_string_to_UTF8(pStr->get(0));
426 Scierror(50, _("%s: Subroutine not found: %s\n"), "odedc", pst);
428 DifferentialEquation::removeDifferentialEquationFunctions();
431 return types::Function::Error;
434 else if (in[iPos]->isList())
436 types::List* pList = in[iPos]->getAs<types::List>();
438 if (pList->getSize() == 0)
440 Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "odedc", iPos + 1, "(string empty)");
441 DifferentialEquation::removeDifferentialEquationFunctions();
444 return types::Function::Error;
447 if (bFuncF && (bFuncJac || pDblNg) && (bFuncG || meth != 3))
449 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "odedc", iPos + 1);
450 DifferentialEquation::removeDifferentialEquationFunctions();
453 return types::Function::Error;
456 if (pList->get(0)->isString())
458 types::String* pStr = pList->get(0)->getAs<types::String>();
464 bOK = deFunctionsManager->setFFunction(pStr);
465 sizeOfpdYData = *YSize;
467 else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
470 bOK = deFunctionsManager->setJacFunction(pStr);
471 if (sizeOfpdYData == 0)
473 sizeOfpdYData = *YSize;
476 else if (bFuncG == false && meth == 3)
479 bOK = deFunctionsManager->setGFunction(pStr);
480 if (sizeOfpdYData == 0)
482 sizeOfpdYData = *YSize;
488 char* pst = wide_string_to_UTF8(pStr->get(0));
489 Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "odedc", iPos + 1, pst);
491 DifferentialEquation::removeDifferentialEquationFunctions();
494 return types::Function::Error;
497 // store in YSize the size of "args" input
498 int* sizeTemp = YSize;
499 int totalSize = sizeOfpdYData;
501 YSize = (int*)MALLOC((sizeOfYSize + pList->getSize() - 1) * sizeof(int));
502 memcpy(YSize, sizeTemp, sizeOfYSize * sizeof(int));
504 std::vector<types::Double*> vpDbl;
505 for (int iter = 0; iter < pList->getSize() - 1; iter++)
507 if (pList->get(iter + 1)->isDouble() == false)
509 Scierror(999, _("%s: Wrong type for input argument #%d: Argument %d in the list must be a matrix.\n"), "odedc", iPos + 1, iter + 1);
512 return types::Function::Error;
515 vpDbl.push_back(pList->get(iter + 1)->getAs<types::Double>());
516 YSize[sizeOfYSize + iter] = vpDbl[iter]->getSize();
517 totalSize += YSize[sizeOfYSize + iter];
520 // store in pdYData "args" input
521 double* pdYDataTemp = pdYData;
522 pdYData = (double*)MALLOC(totalSize * sizeof(double));
523 C2F(dcopy)(&sizeOfpdYData, pdYDataTemp, &one, pdYData, &one);
525 int position = sizeOfpdYData;
526 for (int iter = 0; iter < pList->getSize() - 1; iter++)
528 C2F(dcopy)(&YSize[sizeOfYSize + iter], vpDbl[iter]->get(), &one, &pdYData[position], &one);
529 position += vpDbl[iter]->getSize();
532 sizeOfpdYData = totalSize;
533 sizeOfYSize += pList->getSize() - 1;
537 else if (pList->get(0)->isCallable())
542 deFunctionsManager->setFFunction(pList->get(0)->getAs<types::Callable>());
543 for (int iter = 1; iter < pList->getSize(); iter++)
545 deFunctionsManager->setFArgs(pList->get(iter)->getAs<types::InternalType>());
548 else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
551 deFunctionsManager->setJacFunction(pList->get(0)->getAs<types::Callable>());
552 for (int iter = 1; iter < pList->getSize(); iter++)
554 deFunctionsManager->setJacArgs(pList->get(iter)->getAs<types::InternalType>());
557 else if (bFuncG == false && meth == 3)
560 deFunctionsManager->setGFunction(pList->get(0)->getAs<types::Callable>());
561 for (int iter = 1; iter < pList->getSize(); iter++)
563 deFunctionsManager->setGArgs(pList->get(iter)->getAs<types::InternalType>());
569 Scierror(999, _("%s: Wrong type for input argument #%d: The first argument in the list must be a string or a function.\n"), "odedc", iPos + 1);
570 DifferentialEquation::removeDifferentialEquationFunctions();
573 return types::Function::Error;
578 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a function expected.\n"), "odedc", iPos + 1);
579 DifferentialEquation::removeDifferentialEquationFunctions();
582 return types::Function::Error;
588 int val = (meth == 3) ? 3 : 1;
589 Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "odedc", in.size() + val);
590 DifferentialEquation::removeDifferentialEquationFunctions();
593 return types::Function::Error;
595 if (pDblNg == NULL && meth == 3)
597 Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "odedc", in.size() + 2);
598 DifferentialEquation::removeDifferentialEquationFunctions();
601 return types::Function::Error;
603 if (bFuncG == false && meth == 3)
605 Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "odedc", in.size() + 1);
606 DifferentialEquation::removeDifferentialEquationFunctions();
609 return types::Function::Error;
612 // *** Initialization. ***
616 int istate = 1; // istate = integer flag (input and output). set istate = 1.
617 int iopt = 0; // iopt = 0 to indicate no optional inputs used.
618 int jt = 2; // jacobian type indicator. set jt = 2.
623 double* rwork = NULL;
628 // contain ls0001, lsa001 and eh0001 structures
629 double* dStructTab = NULL;
630 int* iStructTab = NULL;
631 int dStructTabSize = 0;
632 int iStructTabSize = 0;
634 int rwSize = 0; // rwSize = dStructTab + rworkSize
635 int iwSize = 0; // iwSize = iStructTab + iworkSize
637 // structures used by lsoda and lsode
638 double* ls0001d = &(C2F(ls0001).tret);
639 int* ls0001i = &(C2F(ls0001).illin);
640 double* lsa001d = &(C2F(lsa001).tsw);
641 int* lsa001i = &(C2F(lsa001).insufr);
642 double* lsr001d = &(C2F(lsr001).rownr3[0]);
643 int* lsr001i = &(C2F(lsr001).lg0);
644 int* eh0001i = &(C2F(eh0001).mesflg);
647 types::InternalType* pIT = symbol::Context::getInstance()->get(symbol::Symbol(L"%ODEOPTIONS"));
648 if (pIT != NULL && pIT->isDouble())
650 pDblOdeOptions = pIT->getAs<types::Double>();
651 if (pDblOdeOptions->getSize() == 12)
654 itask = (int)pDblOdeOptions->get(0);
655 jt = (int)pDblOdeOptions->get(5);
656 ml = (int)pDblOdeOptions->get(10);
657 mu = (int)pDblOdeOptions->get(11);
659 if (itask != 4 && itask != 5)
670 if (getWarningMode())
672 sciprint(_("%s: Warning: odedc forces itask = %d.\n"), "odedc", itask);
678 if (iopt == 1 && (pDblOdeOptions->get(4) > pDblOdeOptions->get(3))) // hmin > hmax ?
680 Scierror(9999, _("%s: Wrong value of hmin and hmax: hmin = %d is greater than hmax = %d.\n"), "odedc", pDblOdeOptions->get(4), pDblOdeOptions->get(3));
681 DifferentialEquation::removeDifferentialEquationFunctions();
684 return types::Function::Error;
687 if (jt < 0 || jt > 5)
689 Scierror(9999, _("%s: Wrong value of Jacobian type: A number between %d and %d expected.\n"), "odedc", 0, 5);
690 DifferentialEquation::removeDifferentialEquationFunctions();
693 return types::Function::Error;
696 if (iopt == 0 && bFuncJac)
701 if (bFuncJac && (jt == 2 || jt == 5) && getWarningMode())
703 sciprint(_("%s: Warning: Jacobian is given, but not used.\n"), "odedc");
706 if (bFuncJac == false && (jt == 1 || jt == 4))
708 if (getWarningMode())
710 sciprint(_("%s: Warning: No Jacobian external given, but one is required by %ODEOPTIONS(6) value. Jacobian will be estimated.\n"), "odedc");
716 //compute itol and set the tolerances rtol and atol.
722 if (pDblRtol->isScalar())
724 rtol = (double*)MALLOC(sizeof(double));
725 *rtol = pDblRtol->get(0);
729 rtol = pDblRtol->get();
735 rtol = (double*)MALLOC(sizeof(double));
736 if (meth == 6 || meth == 7)
748 if (pDblAtol->isScalar())
750 atol = (double*)MALLOC(sizeof(double));
751 *atol = pDblAtol->get(0);
755 atol = pDblAtol->get();
761 atol = (double*)MALLOC(sizeof(double));
762 if (meth == 6 || meth == 7)
772 // Compute rwork, iwork size and create them.
780 if (pDblW) // structure ls0001 have been restored.
782 nyh = C2F(ls0001).nyh;
787 mxordn = (int)pDblOdeOptions->get(7);
788 mxords = (int)pDblOdeOptions->get(8);
791 if (mxordn > 12 || mxords > 5 || mxordn < 1 || mxords < 1)
793 if (getWarningMode())
795 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);
806 case 3: // lsodar (root)
808 // lrn = 20 + nyh*(mxordn+1) + 3*neq + 3*ng,
809 // lrs = 20 + nyh*(mxords+1) + 3*neq + lmat + 3*ng,
811 // nyh = the initial value of neq,
812 // mxordn = 12, unless a smaller value is given as an
814 // mxords = 5, unless a smaller value is given as an
816 // lmat = length of matrix work space..
817 // lmat = neq**2 + 2 if jt = 1 or 2,
818 // lmat = (2*ml + mu + 1)*neq + 2 if jt = 4 or 5.
820 lrn = 3 * (int)pDblNg->get(0);
822 dStructTabSize = 246 - 241;
823 iStructTabSize = 59 - 50;
824 jroot = (int*)MALLOC((int)pDblNg->get(0) * sizeof(int));
828 // lrn = 20 + nyh*(mxordn+1) + 3*neq,
829 // lrs = 20 + nyh*(mxords+1) + 3*neq + lmat,
831 // nyh = the initial value of neq,
832 // mxordn = 12, unless a smaller value is given as an
834 // mxords = 5, unless a smaller value is given as an
836 // lmat = length of matrix work space..
837 // lmat = neq**2 + 2 if jt = 1 or 2,
838 // lmat = (2*ml + mu + 1)*neq + 2 if jt = 4 or 5.
840 if (jt == 1 || jt == 2)
842 lmat = (sizeYc) * (sizeYc) + 2; // if jt = 1 or 2,
844 else if (jt == 4 || jt == 5)
846 //ml and mu = -1 in all cases
847 lmat = (2 * ml + mu + 1) * (sizeYc) + 2; // if jt = 4 or 5.
850 lrn += 20 + nyh * (mxordn + 1) + 3 * (sizeYc);
851 lrs += 20 + nyh * (mxords + 1) + 3 * (sizeYc) + lmat;
853 rworkSize = max(lrn, lrs);
854 iworkSize = 20 + sizeYc;
856 dStructTabSize += 241;
857 iStructTabSize += 50;
861 case 1: // lsode (adams) (non stiff)
865 case 2: // lsode (stiff)
867 // 20 + nyh*(maxord + 1) + 3*neq + lmat
869 // nyh = the initial value of neq,
870 // maxord = 12 (if meth = 1) or 5 (if meth = 2) (unless a
871 // smaller value is given as an optional input),
872 // lmat = 0 if miter = 0,
873 // lmat = neq**2 + 2 if miter is 1 or 2,
874 // lmat = neq + 2 if miter = 3, and
875 // lmat = (2*ml+mu+1)*neq + 2 if miter is 4 or 5.
877 if (jt == 1 || jt == 2)
879 lmat = (sizeYc) * (sizeYc) + 2;
885 else if (jt == 4 || jt == 5)
887 lmat = (2 * ml + mu + 1) * (sizeYc) + 2;
890 rworkSize = 20 + nyh * (maxord + 1) + 3 * (sizeYc) + lmat;
891 iworkSize = 20; // if jt = 0 or 3
892 if (jt == 1 || jt == 2 || jt == 4 || jt == 5) // iSize = 20 + neq otherwise
894 iworkSize += (sizeYc);
897 dStructTabSize = 219;
901 case 4: // lsdisc (discrete)
907 case 5: // lsrgk (rk)
909 rworkSize = 3 * (sizeYc);
913 case 6: // rkf45 (rkf)
915 rworkSize = 3 + 8 * (sizeYc);
919 case 7: // rksimp (fix)
921 rworkSize = 3 + 8 * (sizeYc);
927 rwSize = rworkSize + dStructTabSize;
928 iwSize = iworkSize + iStructTabSize;
930 rwork = (double*)MALLOC(rworkSize * sizeof(double));
931 iwork = (int*)MALLOC(iworkSize * sizeof(int));
937 if (pDblW->getSize() != rwSize || pDblIw->getSize() != iwSize)
939 Scierror(9999, _("%s: Wrong size for w and iw: w = %d and iw = %d expected.\n"), "odedc", rwSize, iwSize);
940 DifferentialEquation::removeDifferentialEquationFunctions();
947 if (itol == 1 || itol == 3)
955 return types::Function::Error;
958 istate = 2; // 1 means this is the first call | 2 means this is not the first call
960 dStructTab = (double*)MALLOC(dStructTabSize * sizeof(double));
961 iStructTab = (int*)MALLOC(iStructTabSize * sizeof(int));
963 C2F(dcopy)(&rworkSize, pDblW->get(), &one, rwork, &one);
964 C2F(dcopy)(&dStructTabSize, pDblW->get() + rworkSize, &one, dStructTab, &one);
966 for (int i = 0; i < iworkSize; i++)
968 iwork[i] = (int)pDblIw->get(i);
971 for (int i = 0; i < iStructTabSize; i++)
973 iStructTab[i] = (int)pDblIw->get(i + iworkSize);
978 //if iopt = 0 default value are used without set it.
983 rwork[0] = pDblOdeOptions->get(1); // tcrit
986 rwork[4] = pDblOdeOptions->get(2); // h0 =0
987 rwork[5] = pDblOdeOptions->get(3); // hmax = %inf
988 rwork[6] = pDblOdeOptions->get(4); // hmin = 0
989 iwork[0] = (int)pDblOdeOptions->get(10);// ml = -1
990 iwork[1] = (int)pDblOdeOptions->get(11);// mu = -1
991 iwork[4] = (int)pDblOdeOptions->get(9); // ixpr = 0
992 iwork[5] = (int)pDblOdeOptions->get(6); // mxstep = 500
993 iwork[6] = 10; // mxhnil = 10 maximum number of messages printed per problem
994 iwork[7] = (int)pDblOdeOptions->get(7); // mxordn = 12
995 iwork[8] = (int)pDblOdeOptions->get(8); // mxords = 5
1000 if (pDblOdeOptions && pDblOdeOptions->get(9) == 1)
1002 sciprint(_("itask = %d\tmeth = %d\tjactyp = %d\tml = %d\tmu = %d\n"), itask, meth, jt, ml, mu, iopt);
1003 sciprint(_("tcrit = %lf\th0 = %lf\thmax = %lf\thmin = %lf\n"), pDblOdeOptions->get(1), pDblOdeOptions->get(2), pDblOdeOptions->get(3), pDblOdeOptions->get(4));
1006 // get rwork and iwork
1009 if (pDblW && pDblIw)
1015 //restore ls0001 and eh0001 from w (rwork) and iw (iwork).
1016 C2F(dcopy)(&dSize, dStructTab, &one, ls0001d, &one);
1017 memcpy(ls0001i, iStructTab, 39 * sizeof(int));
1022 //restore lsa001 from w (rwork) and iw (iwork).
1023 if (meth == 0 || meth == 3)
1026 C2F(dcopy)(&dSize, dStructTab + dPos, &one, lsa001d, &one);
1027 memcpy(lsa001i, &iStructTab[iPos], 9 * sizeof(int));
1037 C2F(dcopy)(&dSize, dStructTab + dPos, &one, lsr001d, &one);
1038 memcpy(lsr001i, &iStructTab[iPos], 9 * sizeof(int));
1043 // *** Perform operation. ***
1046 double tleft = pDblT0->get(0);
1047 double tright = pDblT0->get(0);
1051 bool bStore = false;
1052 bool bUpdate = false;
1053 bool bOneStep = false;
1055 std::string strMeth;
1058 if (pDblStdel->isScalar() == false)
1060 delta = (int)pDblStdel->get(1);
1063 if (itask == 5 || itask == 3 || itask == 2)
1066 if (getWarningMode() && pDblT->isScalar() == false)
1068 sciprint(_("itask = %d: At most one value of t is allowed, the last element of t is used.\n"), itask);
1074 std::list<double*> pDblYOutList = std::list<double*>();
1075 std::list<double> pDblTOutList = std::list<double>();
1076 int iLastT = pDblT->getSize() - 1;
1077 bool bIntegrateContPart = false;
1081 if (bIntegrateContPart == false)
1083 hf = min(pDblT0->get(0) + (nhpass + delta) * pDblStdel->get(0), pDblT->get(iLastT));
1086 if (fabs(tleft - hf) < 1.0e-12) // update discrete part
1088 bIntegrateContPart = false;
1089 deFunctionsManager->setOdedcFlag();
1091 if (pDblOdeOptions && pDblOdeOptions->get(9) == 1)
1093 sciprint(_("update at t = %lf\n"), tright);
1097 ode_f(&sizeYc, &tright, pdYData, pdYData + sizeYc);
1099 catch (ast::ScilabError &e)
1101 char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
1102 sciprint(_("%s: Update failed at t = %lf\n"), "odedc", tright);
1103 Scierror(999, pstrMsg);
1105 DifferentialEquation::removeDifferentialEquationFunctions();
1120 if (itol == 1 || itol == 3)
1128 return types::Function::Error;
1131 deFunctionsManager->resetOdedcFlag();
1134 double* copy = (double*)MALLOC(*YSize * sizeof(double));
1135 C2F(dcopy)(YSize, pdYData, &one, copy, &one);
1136 pDblYOutList.push_back(copy);
1137 pDblTOutList.push_back(tleft);
1139 else // integrate continuous part
1141 bIntegrateContPart = true;
1143 rwork[0] = hf; // tcrit = hf
1145 if (pDblOdeOptions && pDblOdeOptions->get(9) == 1)
1147 sciprint(_("integ. from tleft= %lf to tf= %lf\n"), tleft, tright);
1154 case 1: // lsode (adams)
1155 case 2: // lsode (stiff)
1158 int jacType = 10 * meth + jt;
1159 C2F(lsode)(ode_f, &sizeYc, pdYData, &tleft, &tright, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jacType);
1165 int ng = (int)pDblNg->get(0);
1166 C2F(lsodar)(ode_f, &sizeYc, pdYData, &tleft, &tright, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt, ode_g, &ng, jroot);
1169 case 4: // lsdisc (discrete)
1172 C2F(lsdisc)(ode_f, &sizeYc, pdYData, &tleft, &tright, rwork, &rworkSize, &istate);
1175 case 5: // lsrgk (rk)
1178 C2F(lsrgk)(ode_f, &sizeYc, pdYData, &tleft, &tright, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1181 case 6: // rkf45 (rkf)
1184 C2F(rkf45)(ode_f, &sizeYc, pdYData, &tleft, &tright, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1187 case 7: // rksimp (fix)
1190 C2F(rksimp)(ode_f, &sizeYc, pdYData, &tleft, &tright, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1193 default: // case 0: // lsoda
1196 C2F(lsoda)(ode_f, &sizeYc, pdYData, &tleft, &tright, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt);
1201 err = checkOdeError(meth, istate);
1202 if (err == 1) // error case
1204 Scierror(999, _("%s: %s exit with state %d.\n"), "odedc", strMeth.c_str(), istate);
1207 catch (ast::ScilabError &e)
1209 char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
1210 sciprint(_("%s: exception caught in '%s' subroutine.\n"), "odedc", strMeth.c_str());
1211 Scierror(999, pstrMsg);
1216 if (err == 1) // error case
1218 DifferentialEquation::removeDifferentialEquationFunctions();
1233 if (itol == 1 || itol == 3)
1241 return types::Function::Error;
1244 if (err == 2) // warning case
1246 if (getWarningMode())
1248 sciprint(_("Integration was stoped at t = %lf.\n"), tleft);
1253 double* copy = (double*)MALLOC(*YSize * sizeof(double));
1254 C2F(dcopy)(YSize, pdYData, &one, copy, &one);
1255 pDblYOutList.push_back(copy);
1256 pDblTOutList.push_back(tleft);
1258 if (meth == 3 && istate == 3 && getWarningMode())
1260 // istate == 3 means the integration was successful, and one or more
1261 // roots were found before satisfying the stop condition
1262 // specified by itask. see jroot.
1264 sciprint(_("%s: Warning: At t = %lf, y is a root, jroot = "), "odedc", tleft);
1265 for (int k = 0; k < pDblNg->get(0); k++)
1267 sciprint("\t%d", jroot[k]);
1274 while (fabs(tleft - pDblT->get(iLastT)) > 1.0e-12 || bIntegrateContPart);
1277 pDblYOut = new types::Double(pDblY0->getRows() + 1, (int)pDblYOutList.size() * pDblY0->getCols());
1278 int iSizeList = (int)pDblYOutList.size();
1279 for (int i = 0; i < iSizeList; i++)
1281 pDblYOut->set(i * (*YSize + 1), pDblTOutList.front());
1282 for (int j = 0; j < *YSize; j++)
1284 pDblYOut->set(i * (*YSize + 1) + (j + 1), pDblYOutList.front()[j]);
1286 pDblYOutList.pop_front();
1287 pDblTOutList.pop_front();
1293 pDblYOut = new types::Double(pDblY0->getRows(), pDblT->getSize() * pDblY0->getCols());
1295 for (int i = 0; i < pDblT->getSize(); i++)
1297 hf = pDblT0->get(0) + (nhpass + delta) * pDblStdel->get(0);
1300 if (pDblOdeOptions && pDblOdeOptions->get(9) == 1)
1302 sciprint("tf - hf = %lf\n", tf - hf);
1305 if (fabs(tf - hf) < 1.0e-12)
1310 if (pDblOdeOptions && pDblOdeOptions->get(9) == 1)
1312 sciprint(_("integ. from tleft= %lf to hf=tf= %lf\n"), tleft, tright);
1320 if (pDblOdeOptions && pDblOdeOptions->get(9) == 1)
1322 sciprint(_("integ. from tleft= %lf to tf= %lf\n"), tleft, tright);
1330 if (pDblOdeOptions && pDblOdeOptions->get(9) == 1)
1332 sciprint(_("integ. from tleft= %lf to hf= %lf\n"), tleft, tright);
1336 rwork[0] = hf; // tcrit = hf
1341 case 1: // lsode (adams)
1342 case 2: // lsode (stiff)
1345 int jacType = 10 * meth + jt;
1346 C2F(lsode)(ode_f, &sizeYc, pdYData, &tleft, &tright, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jacType);
1352 int ng = (int)pDblNg->get(0);
1353 C2F(lsodar)(ode_f, &sizeYc, pdYData, &tleft, &tright, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt, ode_g, &ng, jroot);
1356 case 4: // lsdisc (discrete)
1359 C2F(lsdisc)(ode_f, &sizeYc, pdYData, &tleft, &tright, rwork, &rworkSize, &istate);
1362 case 5: // lsrgk (rk)
1365 C2F(lsrgk)(ode_f, &sizeYc, pdYData, &tleft, &tright, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1368 case 6: // rkf45 (rkf)
1371 C2F(rkf45)(ode_f, &sizeYc, pdYData, &tleft, &tright, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1374 case 7: // rksimp (fix)
1377 C2F(rksimp)(ode_f, &sizeYc, pdYData, &tleft, &tright, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1380 default: // case 0: // lsoda
1383 C2F(lsoda)(ode_f, &sizeYc, pdYData, &tleft, &tright, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt);
1389 err = checkOdeError(meth, istate);
1390 if (err == 1) // error case
1392 Scierror(999, _("%s: %s exit with state %d.\n"), "odedc", strMeth.c_str(), istate);
1395 catch (ast::ScilabError &e)
1397 char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
1398 sciprint(_("%s: exception caught in '%s' subroutine.\n"), "odedc", strMeth.c_str());
1399 Scierror(999, pstrMsg);
1404 if (err == 1) // error case
1406 DifferentialEquation::removeDifferentialEquationFunctions();
1421 if (itol == 1 || itol == 3)
1429 return types::Function::Error;
1434 deFunctionsManager->setOdedcFlag();
1436 if (pDblOdeOptions && pDblOdeOptions->get(9) == 1)
1438 sciprint(_("update at t = %lf\n"), tright);
1443 ode_f(&sizeYc, &tright, pdYData, pdYData + sizeYc);
1445 catch (ast::ScilabError &e)
1447 char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
1448 sciprint(_("%s: Update failed at t = %lf\n"), "odedc", tright);
1449 Scierror(999, pstrMsg);
1450 DifferentialEquation::removeDifferentialEquationFunctions();
1465 if (itol == 1 || itol == 3)
1473 return types::Function::Error;
1476 deFunctionsManager->resetOdedcFlag();
1482 for (int j = 0; j < *YSize; j++)
1484 pDblYOut->set(i * (*YSize) + j, pdYData[j]);
1492 if (err == 2) // warning case
1494 if (getWarningMode())
1496 sciprint(_("Integration was stoped at t = %lf.\n"), tleft);
1499 types::Double* pDblYOutTemp = pDblYOut;
1500 pDblYOut = new types::Double(pDblYOutTemp->getRows(), i + 1);
1501 for (int k = 0; k < pDblYOut->getSize(); k++)
1503 pDblYOut->set(k, pDblYOutTemp->get(k));
1505 delete pDblYOutTemp;
1509 if (meth == 3 && istate == 3 && getWarningMode())
1511 // istate == 3 means the integration was successful, and one or more
1512 // roots were found before satisfying the stop condition
1513 // specified by itask. see jroot.
1515 sciprint(_("%s: Warning: At t = %lf, y is a root, jroot = "), "odedc", tleft);
1516 for (int k = 0; k < pDblNg->get(0); k++)
1518 sciprint("\t%d", jroot[k]);
1522 types::Double* pDblYOutTemp = pDblYOut;
1523 pDblYOut = new types::Double(pDblYOutTemp->getRows(), i + 1);
1524 for (int k = 0; k < pDblYOut->getSize(); k++)
1526 pDblYOut->set(k, pDblYOutTemp->get(k));
1528 delete pDblYOutTemp;
1536 if (_iRetCount > 2) //save ls0001 and eh0001 following w and iw.
1542 if (dStructTab == NULL)
1544 dStructTab = (double*)MALLOC(dStructTabSize * sizeof(double));
1545 iStructTab = (int*)MALLOC(iStructTabSize * sizeof(int));
1548 C2F(dcopy)(&dSize, ls0001d, &one, dStructTab, &one);
1549 memcpy(iStructTab, ls0001i, 39 * sizeof(int));
1555 if (meth == 0 || meth == 3)
1558 C2F(dcopy)(&dSize, lsa001d, &one, dStructTab + dPos, &one);
1559 memcpy(&iStructTab[iPos], lsa001i, 9 * sizeof(int));
1569 C2F(dcopy)(&dSize, lsr001d, &one, dStructTab + dPos, &one);
1570 memcpy(&iStructTab[iPos], lsr001i, 9 * sizeof(int));
1575 memcpy(&iStructTab[iPos], eh0001i, 2 * sizeof(int));
1579 // *** Return result in Scilab. ***
1581 out.push_back(pDblYOut); // y
1583 if (meth == 3 && _iRetCount >= 2)
1588 for (int i = 0; i < pDblNg->get(0); i++)
1596 types::Double* pDblRd = new types::Double(1, sizeOfRd);
1597 //rd: The first entry contains the stopping time.
1598 pDblRd->set(0, C2F(lsr001).tlast);
1599 for (int i = 0; i < pDblNg->get(0); i++)
1604 pDblRd->set(k, (double)i + 1);
1607 out.push_back(pDblRd); // rd
1610 if ((meth < 3 && _iRetCount == 3) || (meth == 3 && _iRetCount == 4))
1612 types::Double* pDblWOut = new types::Double(1, rwSize);
1613 C2F(dcopy)(&rworkSize, rwork, &one, pDblWOut->get(), &one);
1614 C2F(dcopy)(&dStructTabSize, dStructTab, &one, pDblWOut->get() + rworkSize, &one);
1616 types::Double* pDblIwOut = new types::Double(1, iwSize);
1617 for (int i = 0; i < iworkSize; i++)
1619 pDblIwOut->set(i, (double)iwork[i]);
1622 for (int i = 0; i < iStructTabSize; i++)
1624 pDblIwOut->set(iworkSize + i, (double)iStructTab[i]);
1627 out.push_back(pDblWOut); // w
1628 out.push_back(pDblIwOut); // iw
1632 if (itol == 1 || itol == 3) // atol is scalar
1637 if (itol < 3) // rtol is scalar
1658 DifferentialEquation::removeDifferentialEquationFunctions();
1660 return types::Function::OK;