2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2013 - Scilab Enterprises - Paul Bignier
4 * Copyright (C) 2013 - Scilab Enterprises - Cedric DELAMARRE
6 * This file must be used under the terms of the CeCILL.
7 * This source file is licensed as described in the file COPYING, which
8 * you should have received as part of this distribution. The terms
9 * are also available at
10 * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
13 /*--------------------------------------------------------------------------*/
15 #include "differential_equations_gw.hxx"
16 #include "function.hxx"
20 #include "callable.hxx"
21 #include "differentialequationfunctions.hxx"
22 #include "runvisitor.hxx"
23 #include "checkodeerror.hxx"
27 #include "sci_malloc.h"
28 #include "localization.h"
31 #include "scifunctions.h"
32 #include "elem_common.h"
35 /*--------------------------------------------------------------------------*/
36 // This code was inspired by sci_f_daskr.f from the same folder.
37 // daskr() has more optional parameters than daskr(), so the gateway
38 // has to adapt to the new options
39 /*--------------------------------------------------------------------------*/
40 //[y0,nvs[,hotdata]]=daskr(y0,t0,t1[,atol[,rtol]],res[,jac],nh,h[,info[,psol][,pjac]][,hotdata])
41 types::Function::ReturnValue sci_daskr(types::typed_list &in, int _iRetCount, types::typed_list &out)
44 types::Double* pDblX0 = NULL;
45 types::Double* pDblT0 = NULL;
46 types::Double* pDblT = NULL;
47 types::Double* pDblRtol = NULL;
48 types::Double* pDblAtol = NULL;
49 types::Double* pDblHd = NULL;
50 types::Double* pDblNg = NULL;
51 types::Double* pDblE7 = NULL; // 7th element of list info
52 types::Double* pDblE12 = NULL; // 12th element of list info
56 double* pdYData = NULL; // contain y0 following by all args data in list case.
57 double* pdYdotData = NULL;
58 int sizeOfpdYData = 0;
61 int* YSize = NULL; // YSize(1) = size of y0,
62 // YSize(n) = size of Args(n) in list case.
66 int info[20] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
90 // Indicate if the function is given.
92 bool bFuncJac = false;
94 bool bFuncPsol = false;
95 bool bFuncPjac = false;
97 // Indicate if info list is given.
98 bool bListInfo = false;
100 // *** check the minimal number of input args. ***
101 if ((int)in.size() < 6 || (int)in.size() > 13)
103 Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "daskr", 6, 13);
104 return types::Function::Error;
107 // *** check number of output args ***
108 if (_iRetCount != 3 && _iRetCount != 2)
110 Scierror(78, _("%s: Wrong number of output argument(s): %d to %d expected.\n"), "daskr", 2, 3);
111 return types::Function::Error;
114 // *** check type of input args and get it. ***
116 if (in[iPos]->isDouble() == false)
118 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "daskr", iPos + 1);
119 return types::Function::Error;
122 pDblX0 = in[iPos]->getAs<types::Double>();
124 if (pDblX0->isComplex())
126 Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "daskr", iPos + 1);
127 return types::Function::Error;
130 if (pDblX0->getCols() > 2)
132 Scierror(999, _("%s: Wrong size for input argument #%d: A real matrix with %d to %d colomn(s) expected.\n"), "daskr", iPos + 1, 1, 2);
133 return types::Function::Error;
136 if (pDblX0->getCols() == 1)
141 YSize = (int*)MALLOC(sizeOfYSize * sizeof(int));
142 *YSize = pDblX0->getRows();
144 pdYData = (double*)MALLOC(*YSize * sizeof(double));
145 pdYdotData = (double*)MALLOC(*YSize * sizeof(double));
147 C2F(dcopy)(YSize, pDblX0->get(), &iOne, pdYData, &iOne);
148 if (pDblX0->getCols() == 2)
150 C2F(dcopy)(YSize, pDblX0->get() + *YSize, &iOne, pdYdotData, &iOne);
154 memset(pdYdotData, 0x00, *YSize);
159 if (in[iPos]->isDouble() == false)
161 Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "daskr", iPos + 1);
162 return types::Function::Error;
165 pDblT0 = in[iPos]->getAs<types::Double>();
167 if (pDblT0->isScalar() == false)
169 Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), "daskr", iPos + 1);
170 return types::Function::Error;
175 if (in[iPos]->isDouble() == false)
177 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "daskr", iPos + 1);
178 return types::Function::Error;
181 pDblT = in[iPos]->getAs<types::Double>();
183 if (pDblT->isComplex())
185 Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "daskr", iPos + 1);
186 return types::Function::Error;
190 DifferentialEquationFunctions* deFunctionsManager = new DifferentialEquationFunctions(L"daskr");
191 DifferentialEquation::addDifferentialEquationFunctions(deFunctionsManager);
193 deFunctionsManager->setOdeYRows(pDblX0->getRows());
195 for (iPos++; iPos < in.size(); iPos++)
197 if (in[iPos]->isDouble())
199 if (pDblAtol == NULL && bFuncF == false)
201 pDblAtol = in[iPos]->getAs<types::Double>();
202 if (pDblAtol->getSize() != pDblX0->getRows() && pDblAtol->isScalar() == false)
204 Scierror(267, _("%s: Wrong size for input argument #%d: A scalar or a matrix of size %d expected.\n"), "daskr", iPos + 1, pDblX0->getRows());
205 DifferentialEquation::removeDifferentialEquationFunctions();
209 return types::Function::Error;
212 else if (pDblRtol == NULL && bFuncF == false)
214 pDblRtol = in[iPos]->getAs<types::Double>();
215 if (pDblAtol->getSize() != pDblRtol->getSize())
217 Scierror(267, _("%s: Wrong size for input argument #%d: Atol and Rtol must have the same size.\n"), "daskr", iPos + 1, pDblX0->getRows());
218 DifferentialEquation::removeDifferentialEquationFunctions();
222 return types::Function::Error;
225 else if (pDblNg == NULL && bFuncF == true)
227 pDblNg = in[iPos]->getAs<types::Double>();
228 if (pDblNg->isScalar() == false)
230 Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), "daskr", iPos + 1);
231 DifferentialEquation::removeDifferentialEquationFunctions();
235 return types::Function::Error;
237 ng = (int)pDblNg->get(0);
239 else if (pDblHd == NULL && bFuncF == true)
241 pDblHd = in[iPos]->getAs<types::Double>();
242 if ((int)in.size() != iPos + 1)
244 Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "daskr", iPos + 1);
245 DifferentialEquation::removeDifferentialEquationFunctions();
249 return types::Function::Error;
254 Scierror(999, _("%s: Wrong type for input argument #%d: A function expected.\n"), "daskr", iPos + 1);
255 DifferentialEquation::removeDifferentialEquationFunctions();
259 return types::Function::Error;
262 else if (in[iPos]->isCallable())
264 types::Callable* pCall = in[iPos]->getAs<types::Callable>();
267 deFunctionsManager->setFFunction(pCall);
270 else if (bFuncJac == false && pDblNg == NULL)
272 deFunctionsManager->setJacFunction(pCall);
275 else if (bFuncG == false && pDblNg)
277 deFunctionsManager->setGFunction(pCall);
280 else if (bFuncG && bFuncPsol == false)
282 deFunctionsManager->setPsolFunction(pCall);
285 else if (bFuncPsol && bFuncPjac == false)
287 deFunctionsManager->setPjacFunction(pCall);
292 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a list expected.\n"), "daskr", iPos + 1);
293 DifferentialEquation::removeDifferentialEquationFunctions();
297 return types::Function::Error;
300 else if (in[iPos]->isString())
302 types::String* pStr = in[iPos]->getAs<types::String>();
307 bOK = deFunctionsManager->setFFunction(pStr);
310 else if (bFuncJac == false && pDblNg == NULL)
312 bOK = deFunctionsManager->setJacFunction(pStr);
315 else if (bFuncG == false && pDblNg)
317 bOK = deFunctionsManager->setGFunction(pStr);
320 else if (bFuncG && bFuncPsol == false)
322 bOK = deFunctionsManager->setPsolFunction(pStr);
325 else if (bFuncPsol && bFuncPjac == false)
327 bOK = deFunctionsManager->setPjacFunction(pStr);
332 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a list expected.\n"), "daskr", iPos + 1);
333 DifferentialEquation::removeDifferentialEquationFunctions();
337 return types::Function::Error;
342 char* pst = wide_string_to_UTF8(pStr->get(0));
343 Scierror(50, _("%s: Subroutine not found: %s\n"), "daskr", pst);
345 DifferentialEquation::removeDifferentialEquationFunctions();
349 return types::Function::Error;
352 else if (in[iPos]->isList())
354 types::List* pList = in[iPos]->getAs<types::List>();
356 if (pList->getSize() == 0)
358 Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "daskr", iPos + 1, "(string empty)");
359 DifferentialEquation::removeDifferentialEquationFunctions();
363 return types::Function::Error;
366 if (bFuncF && bListInfo)
368 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "daskr", iPos + 1);
369 DifferentialEquation::removeDifferentialEquationFunctions();
373 return types::Function::Error;
376 // functions can be passed in a list
377 if (pList->get(0)->isString())
379 // list(string,... function case
380 types::String* pStr = pList->get(0)->getAs<types::String>();
382 sizeOfpdYData = *YSize;
387 bOK = deFunctionsManager->setFFunction(pStr);
389 else if (bFuncJac == false && pDblNg == NULL)
392 bOK = deFunctionsManager->setJacFunction(pStr);
394 else if (bFuncG == false && pDblNg)
397 bOK = deFunctionsManager->setGFunction(pStr);
399 else if (bFuncG && bFuncPsol == false)
401 bOK = deFunctionsManager->setPsolFunction(pStr);
404 else if (bFuncPsol && bFuncPjac == false)
406 bOK = deFunctionsManager->setPjacFunction(pStr);
412 char* pst = wide_string_to_UTF8(pStr->get(0));
413 Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "daskr", iPos + 1, pst);
415 DifferentialEquation::removeDifferentialEquationFunctions();
419 return types::Function::Error;
422 int* sizeTemp = YSize;
423 int totalSize = sizeOfpdYData;
425 YSize = (int*)MALLOC((sizeOfYSize + pList->getSize() - 1) * sizeof(int));
426 memcpy(YSize, sizeTemp, sizeOfYSize * sizeof(int));
428 std::vector<types::Double*> vpDbl;
429 for (int iter = 0; iter < pList->getSize() - 1; iter++)
431 if (pList->get(iter + 1)->isDouble() == false)
433 Scierror(999, _("%s: Wrong type for input argument #%d: Argument %d in the list must be a matrix.\n"), "daskr", iPos + 1, iter + 1);
434 DifferentialEquation::removeDifferentialEquationFunctions();
438 return types::Function::Error;
441 vpDbl.push_back(pList->get(iter + 1)->getAs<types::Double>());
442 YSize[sizeOfYSize + iter] = vpDbl[iter]->getSize();
443 totalSize += YSize[sizeOfYSize + iter];
446 double* pdYDataTemp = pdYData;
447 pdYData = (double*)MALLOC(totalSize * sizeof(double));
448 C2F(dcopy)(&sizeOfpdYData, pdYDataTemp, &iOne, pdYData, &iOne);
450 int position = sizeOfpdYData;
451 for (int iter = 0; iter < pList->getSize() - 1; iter++)
453 C2F(dcopy)(&YSize[sizeOfYSize + iter], vpDbl[iter]->get(), &iOne, &pdYData[position], &iOne);
454 position += vpDbl[iter]->getSize();
457 sizeOfpdYData = totalSize;
458 sizeOfYSize += pList->getSize() - 1;
462 else if (pList->get(0)->isCallable())
464 // list(macro,... function case
468 deFunctionsManager->setFFunction(pList->get(0)->getAs<types::Callable>());
469 for (int iter = 1; iter < pList->getSize(); iter++)
471 deFunctionsManager->setFArgs(pList->get(iter)->getAs<types::InternalType>());
474 else if (bFuncJac == false && pDblNg == NULL)
477 deFunctionsManager->setJacFunction(pList->get(0)->getAs<types::Callable>());
478 for (int iter = 1; iter < pList->getSize(); iter++)
480 deFunctionsManager->setJacArgs(pList->get(iter)->getAs<types::InternalType>());
483 else if (bFuncG == false && pDblNg)
486 deFunctionsManager->setGFunction(pList->get(0)->getAs<types::Callable>());
487 for (int iter = 1; iter < pList->getSize(); iter++)
489 deFunctionsManager->setGArgs(pList->get(iter)->getAs<types::InternalType>());
492 else if (bFuncG && bFuncPsol == false)
495 deFunctionsManager->setPsolFunction(pList->get(0)->getAs<types::Callable>());
496 for (int iter = 1; iter < pList->getSize(); iter++)
498 deFunctionsManager->setPsolArgs(pList->get(iter)->getAs<types::InternalType>());
501 else if (bFuncPsol && bFuncPjac == false)
504 deFunctionsManager->setPjacFunction(pList->get(0)->getAs<types::Callable>());
505 for (int iter = 1; iter < pList->getSize(); iter++)
507 deFunctionsManager->setPjacArgs(pList->get(iter)->getAs<types::InternalType>());
511 else if (pList->get(0)->isDouble() && bFuncF == true)
513 // list(double,... then this list is the info argument
514 if (pList->getSize() != 14)
516 Scierror(267, _("%s: Wrong size for input argument #%d: A list of size %d expected.\n"), "daskr", iPos + 1, 7);
517 DifferentialEquation::removeDifferentialEquationFunctions();
521 return types::Function::Error;
524 for (int i = 0; i < 14; i++) // info = list([],0,[],[],[],0,[],0,[],0,[],[],[],1)
526 if (pList->get(i)->isDouble() == false || (pList->get(i)->getAs<types::Double>()->isScalar() == false && (i == 1 || i == 5 || i == 7 || i == 9 || i == 13)))
528 if (i == 1 || i == 5 || i == 7 || i == 9 || i == 13)
530 Scierror(999, _("%s: Wrong type for input argument #%d: Element %d in the info list must be a scalar.\n"), "daskr", iPos + 1, i);
534 Scierror(999, _("%s: Wrong type for input argument #%d: Element %d in the info list must be a matrix.\n"), "daskr", iPos + 1, i);
536 DifferentialEquation::removeDifferentialEquationFunctions();
540 return types::Function::Error;
544 // -- subvariable tstop(info) --
545 types::Double* pDblTemp = pList->get(0)->getAs<types::Double>();
546 if (pDblTemp->getSize() != 0)
549 tstop = pDblTemp->get(0);
552 // -- subvariable imode(info) --
553 info[2] = (int)pList->get(1)->getAs<types::Double>()->get(0);
555 // -- subvariable band(info) --
556 pDblTemp = pList->get(2)->getAs<types::Double>();
557 if (pDblTemp->getSize() == 2)
560 ml = (int)pDblTemp->get(0);
561 mu = (int)pDblTemp->get(1);
562 deFunctionsManager->setMl(ml);
563 deFunctionsManager->setMu(mu);
565 else if (pDblTemp->getSize() != 0)
567 Scierror(267, _("%s: Wrong size for input argument #%d: Argument %d in te list must be of size %d.\n"), "daskr", iPos + 1, 3, 2);
568 DifferentialEquation::removeDifferentialEquationFunctions();
572 return types::Function::Error;
575 // -- subvariable maxstep(info) --
576 pDblTemp = pList->get(3)->getAs<types::Double>();
577 if (pDblTemp->getSize() != 0)
580 maxstep = pDblTemp->get(0);
583 // -- subvariable stepin(info) --
584 pDblTemp = pList->get(4)->getAs<types::Double>();
585 if (pDblTemp->getSize() != 0)
588 stepin = pDblTemp->get(0);
591 // -- subvariable nonneg(info) --
592 info[9] = (int)pList->get(5)->getAs<types::Double>()->get(0);
594 // -- subvariable consistent(info) --
595 pDblTemp = pList->get(6)->getAs<types::Double>();
596 if (pDblTemp->getSize() != 0)
599 if (info[9] == 0 || info[9] == 2)
605 LID = 40 + pDblX0->getRows();
610 // -- subvariable iteration(info) --
611 if (pList->get(7)->getAs<types::Double>()->get(0) == 1)
616 // -- subvariable defaultKrylov(info) --
617 pDblTemp = pList->get(8)->getAs<types::Double>();
618 if (pDblTemp->getSize() == 0)
620 // maxl and kmp need default values
621 maxl = min(5, pDblX0->getRows());
626 // info then looks like list(..., [maxl kmp nrmax epli],...)
628 maxl = (int)pDblTemp->get(0);
629 kmp = (int)pDblTemp->get(1);
630 nrmax = (int)pDblTemp->get(2);
631 epli = pDblTemp->get(3);
634 // -- subvariable justConsistentComp(info) --
635 dTemp = pList->get(9)->getAs<types::Double>()->get(0);
638 // Check that info(11) = 1, meaning that the provided initial values
639 // are not consistent
646 // -- subvariable psolJac(info) --
647 dTemp = pList->get(10)->getAs<types::Double>()->get(0);
653 // -- subvariable excludeAlgebraic(info) --
654 pDblTemp = pList->get(11)->getAs<types::Double>();
655 if (pDblTemp->getSize() != 0)
658 if (info[9] == 0 || info[9] == 2)
664 LID = 40 + pDblX0->getRows();
669 // -- subvariable defaultHeuristic(info) --
670 pDblTemp = pList->get(12)->getAs<types::Double>();
671 if (pDblTemp->getSize() != 0)
673 // info then looks like list(..., [mxnit mxnj lsoff stptol epinit],...)
675 mxnit = (int)pDblTemp->get(0);
676 mxnj = (int)pDblTemp->get(1);
677 mxnh = (int)pDblTemp->get(2);
678 lsoff = (int)pDblTemp->get(3);
679 stptol = pDblTemp->get(4);
680 epinit = pDblTemp->get(5);
683 // -- subvariable verbosity(info) --
684 dTemp = pList->get(13)->getAs<types::Double>()->get(0);
685 if (dTemp == 1 || dTemp == 2)
687 info[17] = (int)dTemp;
694 Scierror(999, _("%s: Wrong type for input argument #%d: The first argument in the list must be a string, a function or a matrix in case of argument info.\n"), "daskr", iPos + 1);
695 DifferentialEquation::removeDifferentialEquationFunctions();
699 return types::Function::Error;
704 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a function expected.\n"), "daskr", iPos + 1);
705 DifferentialEquation::removeDifferentialEquationFunctions();
709 return types::Function::Error;
715 Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "daskr", (int)in.size() + 3);
716 DifferentialEquation::removeDifferentialEquationFunctions();
720 return types::Function::Error;
725 Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "daskr", (int)in.size() + 2);
726 DifferentialEquation::removeDifferentialEquationFunctions();
730 return types::Function::Error;
735 Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "daskr", (int)in.size() + 1);
736 DifferentialEquation::removeDifferentialEquationFunctions();
740 return types::Function::Error;
743 if (bFuncJac == true)
748 // *** Initialization. ***
749 double t0 = pDblT0->get(0);
753 //compute itol and set the tolerances rtol and atol.
756 double rpar[2] = {0, 0 };
758 memset(ipar, 0x00, 30 * sizeof(int));
762 if (pDblAtol->isScalar())
764 atol = (double*)MALLOC(sizeof(double));
765 *atol = pDblAtol->get(0);
769 atol = pDblAtol->get();
775 atol = (double*)MALLOC(sizeof(double));
781 if (pDblRtol->isScalar())
783 rtol = (double*)MALLOC(sizeof(double));
784 *rtol = pDblRtol->get(0);
788 rtol = pDblRtol->get();
791 else // if rtol is not given atol will be used as a scalar.
793 if (pDblAtol && pDblAtol->isScalar() == false) // info[1] == 1
795 double dblSrc = 1.e-9;
796 int iSize = pDblAtol->getSize();
799 rtol = (double*)MALLOC(iSize * sizeof(double));
800 C2F(dcopy)(&iSize, &dblSrc, &iZero, rtol, &iOne);
804 rtol = (double*)MALLOC(sizeof(double));
809 // Compute rwork, iwork size.
813 double* rwork = NULL;
816 // compute size of rwork
820 rworksize += max(maxord + 4, 7) * pDblX0->getRows() + 3 * ng;
823 // For the full (dense) JACOBIAN case
824 rworksize += pDblX0->getRows() * pDblX0->getRows();
826 else if (info[4] == 1)
828 // For the banded user-defined JACOBIAN case
829 rworksize += (2 * ml + mu + 1) * pDblX0->getRows();
831 else if (info[4] == 0)
833 // For the banded finite-difference-generated JACOBIAN case
834 rworksize += (2 * ml + mu + 1) * pDblX0->getRows() + 2 * (pDblX0->getRows() / (ml + mu + 1) + 1);
837 else // info[11] == 1
839 // LENWP is the length ot the rwork segment containing
840 // the matrix elements of the preconditioner P
842 LENWP = pDblX0->getRows() * pDblX0->getRows();
843 rworksize += (maxord + 5) * pDblX0->getRows() + 3 * ng
844 + (maxl + 3 + min(1, maxl - kmp)) * pDblX0->getRows()
845 + (maxl + 3) * maxl + 1 + LENWP;
850 rworksize += pDblX0->getRows();
853 // compute size of iwork
854 // liw = 40, then augment size according to the case
858 iworksize += pDblX0->getRows();
860 else // info[11] == 1
862 // LENIWP is the length ot the iwork segment containing
863 // the matrix indexes of the preconditioner P
864 // (compressed sparse row format)
865 LENIWP = 2 * pDblX0->getRows() * pDblX0->getRows();
869 if (info[9] == 1 || info[9] == 3)
871 iworksize += pDblX0->getRows();
874 if (info[10] == 1 || info[15] == 1)
876 iworksize += pDblX0->getRows();
879 iwork = (int*)MALLOC(iworksize * sizeof(int));
880 rwork = (double*)MALLOC(rworksize * sizeof(double));
881 root = (int*)MALLOC(ng * sizeof(int));
885 if (iworksize + rworksize != pDblHd->getSize())
887 Scierror(77, _("%s: Wrong size for input argument(s) %d: %d expected.\n"), "daskr", in.size(), iworksize + rworksize);
888 DifferentialEquation::removeDifferentialEquationFunctions();
894 if (pDblAtol == NULL || pDblAtol->isScalar())
898 if (pDblRtol == NULL || pDblRtol->isScalar())
902 return types::Function::Error;
905 C2F(dcopy)(&rworksize, pDblHd->get(), &iOne, rwork, &iOne);
907 for (int i = 0; i < iworksize; i++)
909 iwork[i] = (int)pDblHd->get(rworksize + i);
938 for (int i = 0; i < pDblX0->getRows(); i++)
940 iwork[LID + i] = pDblE7->get(i);
961 ipar[0] = pDblX0->getRows();
962 ipar[1] = pDblX0->getRows();
972 for (int i = 0; i < pDblX0->getRows(); i++)
974 iwork[LID + i] = pDblE12->get(i);
988 iwork[16] = rworksize;
989 iwork[17] = iworksize;
991 // *** Perform operation. ***
992 std::list<types::Double*> lpDblOut;
993 int size = pDblX0->getRows();
994 int rowsOut = 1 + pDblX0->getRows() * 2;
998 for (int i = 0; i < pDblT->getSize(); i++)
1000 types::Double* pDblOut = new types::Double(rowsOut, 1);
1001 lpDblOut.push_back(pDblOut);
1003 double t = pDblT->get(i);
1005 pDblOut->set(pos, t);
1010 C2F(dcopy)(&size, pdYData, &iOne, pDblOut->get() + pos, &iOne);
1011 pos += pDblX0->getRows();
1012 C2F(dcopy)(&size, pdYdotData, &iOne, pDblOut->get() + pos, &iOne);
1019 // SUBROUTINE DDASKR(RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL,
1020 // IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC, PSOL,
1022 C2F(ddaskr)(dassl_f, YSize, &t0, pdYData, pdYdotData, &t,
1023 info, rtol, atol, &idid, rwork, &rworksize,
1024 iwork, &iworksize, rpar, ipar,
1025 info[14] == 1 ? (void*)daskr_pjac : (void*)dassl_jac,
1026 daskr_psol, dasrt_g, &ng, root);
1028 // values of idid says the same things that in ddasrt function,
1029 // except these two values.
1031 // The integration to TSTOP was successfully completed (T=TSTOP)
1037 // One or more root found (eq to idid = 4 for dassl)
1048 // use the same error function that dasrt
1049 iret = checkError(idid, "daskr");
1050 if (iret == 1) // error
1052 Scierror(999, _("%s: ddaskr return with state %d.\n"), "daskr", ididtmp);
1055 catch (ast::ScilabError &e)
1057 char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
1058 sciprint(_("%s: exception caught in '%s' subroutine.\n"), "daskr", "ddaskr");
1059 Scierror(999, pstrMsg);
1060 // set iret to 1 for free allocated memory
1064 // free memory when error occurred
1068 DifferentialEquation::removeDifferentialEquationFunctions();
1075 if (pDblAtol == NULL || pDblAtol->isScalar())
1079 if (pDblRtol == NULL || pDblRtol->isScalar())
1083 return types::Function::Error;
1087 C2F(dcopy)(&size, pdYData, &iOne, pDblOut->get() + pos, &iOne);
1089 C2F(dcopy)(&size, pdYdotData, &iOne, pDblOut->get() + pos, &iOne);
1091 if (iret == 2) // warning
1093 if (ididtmp == -1 || ididtmp == 5 || ididtmp == 6)
1095 pDblOut->set(0, t0);
1100 // iret == 0, check idid flag
1103 pDblOut->set(0, t0);
1106 else if (idid == -2)
1119 // *** Return result in Scilab. ***
1120 types::Double* pDblOut = new types::Double(rowsOut, (int)lpDblOut.size());
1122 int sizeOfList = (int)lpDblOut.size();
1123 for (int i = 0; i < sizeOfList; i++)
1125 int pos = i * rowsOut;
1126 C2F(dcopy)(&rowsOut, lpDblOut.front()->get(), &iOne, pDblOut->get() + pos, &iOne);
1127 lpDblOut.pop_front();
1129 out.push_back(pDblOut);
1132 for (int i = 0; i < ng; i++)
1139 types::Double* pDblRoot = new types::Double(1, sizeOfRoot);
1140 pDblRoot->set(0, t0);
1142 for (int i = 0; i < ng; i++)
1147 pDblRoot->set(j, i + 1);
1150 out.push_back(pDblRoot);
1152 if (_iRetCount == 3)
1154 types::Double* pDblHdOut = new types::Double(rworksize + iworksize, 1);
1155 C2F(dcopy)(&rworksize, rwork, &iOne, pDblHdOut->get(), &iOne);
1157 for (int i = 0; i < iworksize; i++)
1159 pDblHdOut->set(rworksize + i, (double)iwork[i]);
1162 out.push_back(pDblHdOut);
1165 // *** FREE memory. ***
1166 if (pDblAtol == NULL || pDblAtol->isScalar())
1171 if (pDblRtol == NULL || pDblRtol->isScalar())
1183 DifferentialEquation::removeDifferentialEquationFunctions();
1185 return types::Function::OK;