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 "checkodeerror.hxx"
26 #include "sci_malloc.h"
27 #include "localization.h"
30 #include "scifunctions.h"
31 #include "elem_common.h"
35 /*--------------------------------------------------------------------------*/
36 types::Function::ReturnValue sci_dassl(types::typed_list &in, int _iRetCount, types::typed_list &out)
39 types::Double* pDblX0 = NULL;
40 types::Double* pDblT0 = NULL;
41 types::Double* pDblT = NULL;
42 types::Double* pDblRtol = NULL;
43 types::Double* pDblAtol = NULL;
44 types::Double* pDblHd = NULL;
47 double* pdYData = NULL; // contain y0 following by all args data in list case.
48 double* pdYdotData = NULL;
49 int sizeOfpdYData = 0;
54 int info[15] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
63 int* YSize = NULL; // YSize(1) = size of y0,
64 // YSize(n) = size of Args(n) in list case.
66 // Indicate if the function is given.
68 bool bFuncJac = false;
70 // Indicate if info list is given.
71 bool bListInfo = false;
73 // *** check the minimal number of input args. ***
74 if (in.size() < 4 || in.size() > 9)
76 Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "dassl", 4, 9);
77 return types::Function::Error;
80 // *** check number of output args ***
83 Scierror(78, _("%s: Wrong number of output argument(s): %d to %d expected.\n"), "dassl", 1, 2);
84 return types::Function::Error;
87 // *** check type of input args and get it. ***
89 if (in[iPos]->isDouble() == false)
91 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "dassl", iPos + 1);
92 return types::Function::Error;
95 pDblX0 = in[iPos]->getAs<types::Double>();
97 if (pDblX0->isComplex())
99 Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "dassl", iPos + 1);
100 return types::Function::Error;
103 if (pDblX0->getCols() > 2)
105 Scierror(999, _("%s: Wrong size for input argument #%d: A real matrix with %d to %d colomn(s) expected.\n"), "dassl", iPos + 1, 1, 2);
106 return types::Function::Error;
109 if (pDblX0->getCols() == 1)
116 if (in[iPos]->isDouble() == false)
118 Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "dassl", iPos + 1);
119 return types::Function::Error;
122 pDblT0 = in[iPos]->getAs<types::Double>();
124 if (pDblT0->isScalar() == false)
126 Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), "dassl", iPos + 1);
127 return types::Function::Error;
132 if (in[iPos]->isDouble() == false)
134 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "dassl", iPos + 1);
135 return types::Function::Error;
138 pDblT = in[iPos]->getAs<types::Double>();
140 if (pDblT->isComplex())
142 Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "dassl", iPos + 1);
143 return types::Function::Error;
147 DifferentialEquationFunctions* deFunctionsManager = new DifferentialEquationFunctions(L"dassl");
148 DifferentialEquation::addDifferentialEquationFunctions(deFunctionsManager);
150 YSize = (int*)MALLOC(sizeOfYSize * sizeof(int));
151 *YSize = pDblX0->getRows();
153 pdYData = (double*)MALLOC(*YSize * sizeof(double));
154 pdYdotData = (double*)MALLOC(*YSize * sizeof(double));
156 C2F(dcopy)(YSize, pDblX0->get(), &one, pdYData, &one);
157 if (pDblX0->getCols() == 2)
159 C2F(dcopy)(YSize, pDblX0->get() + *YSize, &one, pdYdotData, &one);
163 memset(pdYdotData, 0x00, *YSize);
166 deFunctionsManager->setOdeYRows(pDblX0->getRows());
168 for (iPos++; iPos < in.size(); iPos++)
170 if (in[iPos]->isDouble())
172 if (pDblAtol == NULL && bFuncF == false)
174 pDblAtol = in[iPos]->getAs<types::Double>();
175 if (pDblAtol->getSize() != pDblX0->getRows() && pDblAtol->isScalar() == false)
177 Scierror(267, _("%s: Wrong size for input argument #%d: A scalar or a matrix of size %d expected.\n"), "dassl", iPos + 1, pDblX0->getRows());
178 DifferentialEquation::removeDifferentialEquationFunctions();
182 return types::Function::Error;
185 else if (pDblRtol == NULL && bFuncF == false)
187 pDblRtol = in[iPos]->getAs<types::Double>();
188 if (pDblAtol->getSize() != pDblRtol->getSize())
190 Scierror(267, _("%s: Wrong size for input argument #%d: Atol and Rtol must have the same size.\n"), "dassl", iPos + 1, pDblX0->getRows());
191 DifferentialEquation::removeDifferentialEquationFunctions();
195 return types::Function::Error;
198 else if (pDblHd == NULL && bFuncF == true)
200 pDblHd = in[iPos]->getAs<types::Double>();
201 if (in.size() != iPos + 1)
203 Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "dassl", iPos + 1);
204 DifferentialEquation::removeDifferentialEquationFunctions();
208 return types::Function::Error;
213 Scierror(999, _("%s: Wrong type for input argument #%d: A function expected.\n"), "dassl", iPos + 1);
214 DifferentialEquation::removeDifferentialEquationFunctions();
218 return types::Function::Error;
221 else if (in[iPos]->isCallable())
223 types::Callable* pCall = in[iPos]->getAs<types::Callable>();
226 deFunctionsManager->setFFunction(pCall);
229 else if (bFuncJac == false)
231 deFunctionsManager->setJacFunction(pCall);
236 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a list expected.\n"), "dassl", iPos + 1);
237 DifferentialEquation::removeDifferentialEquationFunctions();
241 return types::Function::Error;
244 else if (in[iPos]->isString())
246 types::String* pStr = in[iPos]->getAs<types::String>();
251 bOK = deFunctionsManager->setFFunction(pStr);
254 else if (bFuncJac == false)
256 bOK = deFunctionsManager->setJacFunction(pStr);
261 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a list expected.\n"), "dassl", iPos + 1);
262 DifferentialEquation::removeDifferentialEquationFunctions();
266 return types::Function::Error;
271 char* pst = wide_string_to_UTF8(pStr->get(0));
272 Scierror(50, _("%s: Subroutine not found: %s\n"), "dassl", pst);
274 DifferentialEquation::removeDifferentialEquationFunctions();
278 return types::Function::Error;
281 else if (in[iPos]->isList())
283 types::List* pList = in[iPos]->getAs<types::List>();
285 if (pList->getSize() == 0)
287 Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "dassl", iPos + 1, "(string empty)");
288 DifferentialEquation::removeDifferentialEquationFunctions();
292 return types::Function::Error;
295 if (bFuncF && bListInfo)
297 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "dassl", iPos + 1);
298 DifferentialEquation::removeDifferentialEquationFunctions();
302 return types::Function::Error;
305 if (pList->get(0)->isString())
307 types::String* pStr = pList->get(0)->getAs<types::String>();
313 bOK = deFunctionsManager->setFFunction(pStr);
314 sizeOfpdYData = *YSize;
316 else if (bFuncJac == false)
319 bOK = deFunctionsManager->setJacFunction(pStr);
320 if (sizeOfpdYData == 0)
322 sizeOfpdYData = *YSize;
328 char* pst = wide_string_to_UTF8(pStr->get(0));
329 Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "dassl", iPos + 1, pst);
331 DifferentialEquation::removeDifferentialEquationFunctions();
335 return types::Function::Error;
338 int* sizeTemp = YSize;
339 int totalSize = sizeOfpdYData;
341 YSize = (int*)MALLOC((sizeOfYSize + pList->getSize() - 1) * sizeof(int));
342 memcpy(YSize, sizeTemp, sizeOfYSize * sizeof(int));
344 std::vector<types::Double*> vpDbl;
345 for (int iter = 0; iter < pList->getSize() - 1; iter++)
347 if (pList->get(iter + 1)->isDouble() == false)
349 Scierror(999, _("%s: Wrong type for input argument #%d: Argument %d in the list must be a matrix.\n"), "dassl", iPos + 1, iter + 1);
350 DifferentialEquation::removeDifferentialEquationFunctions();
354 return types::Function::Error;
357 vpDbl.push_back(pList->get(iter + 1)->getAs<types::Double>());
358 YSize[sizeOfYSize + iter] = vpDbl[iter]->getSize();
359 totalSize += YSize[sizeOfYSize + iter];
362 double* pdYDataTemp = pdYData;
363 pdYData = (double*)MALLOC(totalSize * sizeof(double));
364 C2F(dcopy)(&sizeOfpdYData, pdYDataTemp, &one, pdYData, &one);
366 int position = sizeOfpdYData;
367 for (int iter = 0; iter < pList->getSize() - 1; iter++)
369 C2F(dcopy)(&YSize[sizeOfYSize + iter], vpDbl[iter]->get(), &one, &pdYData[position], &one);
370 position += vpDbl[iter]->getSize();
373 sizeOfpdYData = totalSize;
374 sizeOfYSize += pList->getSize() - 1;
378 else if (pList->get(0)->isCallable())
383 deFunctionsManager->setFFunction(pList->get(0)->getAs<types::Callable>());
384 for (int iter = 1; iter < pList->getSize(); iter++)
386 deFunctionsManager->setFArgs(pList->get(iter)->getAs<types::InternalType>());
389 else if (bFuncJac == false)
392 deFunctionsManager->setJacFunction(pList->get(0)->getAs<types::Callable>());
393 for (int iter = 1; iter < pList->getSize(); iter++)
395 deFunctionsManager->setJacArgs(pList->get(iter)->getAs<types::InternalType>());
399 else if (pList->get(0)->isDouble() && bFuncF == true)
401 if (pList->getSize() != 7)
403 Scierror(267, _("%s: Wrong size for input argument #%d: A list of size %d expected.\n"), "dassl", iPos + 1, 7);
404 DifferentialEquation::removeDifferentialEquationFunctions();
408 return types::Function::Error;
411 for (int i = 0; i < 7; i++) // info = list([],0,[],[],[],0,0)
413 if (pList->get(i)->isDouble() == false || (pList->get(i)->getAs<types::Double>()->isScalar() == false && (i == 1 || i == 5 || i == 6)))
415 if (i == 1 || i == 5 || i == 6)
417 Scierror(999, _("%s: Wrong type for input argument #%d: Element %d in the info list must be a scalar.\n"), "dassl", iPos + 1, i);
421 Scierror(999, _("%s: Wrong type for input argument #%d: Element %d in the info list must be a matrix.\n"), "dassl", iPos + 1, i);
423 DifferentialEquation::removeDifferentialEquationFunctions();
427 return types::Function::Error;
431 types::Double* pDblTemp = pList->get(0)->getAs<types::Double>();
432 if (pDblTemp->getSize() != 0)
435 tstop = pDblTemp->get(0);
438 info[2] = (int)pList->get(1)->getAs<types::Double>()->get(0);
440 pDblTemp = pList->get(2)->getAs<types::Double>();
441 if (pDblTemp->getSize() == 2)
444 ml = (int)pDblTemp->get(0);
445 mu = (int)pDblTemp->get(1);
446 deFunctionsManager->setMl(ml);
447 deFunctionsManager->setMu(mu);
449 else if (pDblTemp->getSize() != 0)
451 Scierror(267, _("%s: Wrong size for input argument #%d: Argument %d in te list must be of size %d.\n"), "dassl", iPos + 1, 3, 2);
452 DifferentialEquation::removeDifferentialEquationFunctions();
456 return types::Function::Error;
459 pDblTemp = pList->get(3)->getAs<types::Double>();
460 if (pDblTemp->getSize() != 0)
463 maxstep = pDblTemp->get(0);
466 pDblTemp = pList->get(4)->getAs<types::Double>();
467 if (pDblTemp->getSize() != 0)
470 stepin = pDblTemp->get(0);
473 info[9] = (int)pList->get(5)->getAs<types::Double>()->get(0);
474 if (pList->get(6)->getAs<types::Double>()->get(0) == 1)
483 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"), "dassl", iPos + 1);
484 DifferentialEquation::removeDifferentialEquationFunctions();
488 return types::Function::Error;
493 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a function expected.\n"), "dassl", iPos + 1);
494 DifferentialEquation::removeDifferentialEquationFunctions();
498 return types::Function::Error;
504 Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "dassl", in.size() + 1);
505 DifferentialEquation::removeDifferentialEquationFunctions();
509 return types::Function::Error;
512 if (bFuncJac == true)
517 // *** Initialization. ***
518 double t0 = pDblT0->get(0);
524 //compute itol and set the tolerances rtol and atol.
530 if (pDblAtol->isScalar())
532 atol = (double*)MALLOC(sizeof(double));
533 *atol = pDblAtol->get(0);
537 atol = pDblAtol->get();
543 atol = (double*)MALLOC(sizeof(double));
549 if (pDblRtol->isScalar())
551 rtol = (double*)MALLOC(sizeof(double));
552 *rtol = pDblRtol->get(0);
556 rtol = pDblRtol->get();
559 else // if rtol is not given atol will be used as a scalar.
561 if (pDblAtol && pDblAtol->isScalar() == false) // info[1] == 1
563 double dblSrc = 1.e-9;
564 int iSize = pDblAtol->getSize();
568 rtol = (double*)MALLOC(iSize * sizeof(double));
569 C2F(dcopy)(&iSize, &dblSrc, &iZero, rtol, &iOne);
573 rtol = (double*)MALLOC(sizeof(double));
578 // Compute rwork, iwork size.
580 int iworksize = 20 + pDblX0->getRows();
583 double* rwork = NULL;
587 rworksize = 40 + (maxord + 4) * pDblX0->getRows() + pDblX0->getRows() * pDblX0->getRows();
589 else if (info[4] == 1)
591 rworksize = 40 + (maxord + 4) * pDblX0->getRows() + (2 * ml + mu + 1) * pDblX0->getRows();
593 else if (info[4] == 0)
595 rworksize = 40 + (maxord + 4) * pDblX0->getRows() + (2 * ml + mu + 1) * pDblX0->getRows() + 2 * (pDblX0->getRows() / (ml + mu + 1) + 1);
598 iwork = (int*)MALLOC(iworksize * sizeof(int));
599 rwork = (double*)MALLOC(rworksize * sizeof(double));
603 if (iworksize + rworksize != pDblHd->getSize())
605 Scierror(77, _("%s: Wrong size for input argument(s) %d: %d expected.\n"), "dassl", in.size(), iworksize + rworksize);
606 DifferentialEquation::removeDifferentialEquationFunctions();
612 if (pDblAtol == NULL || pDblAtol->isScalar())
616 if (pDblRtol == NULL || pDblRtol->isScalar())
620 return types::Function::Error;
623 C2F(dcopy)(&rworksize, pDblHd->get(), &one, rwork, &one);
625 for (int i = 0; i < iworksize; i++)
627 iwork[i] = (int)pDblHd->get(rworksize + i);
654 // *** Perform operation. ***
655 std::list<types::Double*> lpDblOut;
656 int size = pDblX0->getRows();
657 int rowsOut = 1 + pDblX0->getRows() * 2;
660 for (int i = 0; i < pDblT->getSize(); i++)
662 types::Double* pDblOut = new types::Double(rowsOut, 1);
663 lpDblOut.push_back(pDblOut);
665 double t = pDblT->get(i);
667 pDblOut->set(pos, t);
672 C2F(dcopy)(&size, pdYData, &one, pDblOut->get() + pos, &one);
673 pos += pDblX0->getRows();
674 C2F(dcopy)(&size, pdYdotData, &one, pDblOut->get() + pos, &one);
681 C2F(dassl)(dassl_f, YSize, &t0, pdYData, pdYdotData, &t, info, rtol, atol, &idid, rwork, &rworksize, iwork, &iworksize, &rpar, &ipar, dassl_jac);
683 iret = checkError(idid, "dassl");
684 if (iret == 1) // error
686 Scierror(999, _("%s: dassl return with state %d.\n"), "dassl", idid);
689 catch (ast::ScilabError &e)
691 char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
692 sciprint(_("%s: exception caught in '%s' subroutine.\n"), "dassl", "dassl");
693 Scierror(999, pstrMsg);
694 // set iret to 1 for FREE allocated memory
698 if (iret == 1) // error
701 DifferentialEquation::removeDifferentialEquationFunctions();
707 if (pDblAtol == NULL || pDblAtol->isScalar())
711 if (pDblRtol == NULL || pDblRtol->isScalar())
715 return types::Function::Error;
719 C2F(dcopy)(&size, pdYData, &one, pDblOut->get() + pos, &one);
721 C2F(dcopy)(&size, pdYdotData, &one, pDblOut->get() + pos, &one);
723 if (iret == 2) // warning
747 // *** Return result in Scilab. ***
748 types::Double* pDblOut = new types::Double(rowsOut, (int)lpDblOut.size());
750 int sizeOfList = (int)lpDblOut.size();
751 for (int i = 0; i < sizeOfList; i++)
753 int pos = i * rowsOut;
754 C2F(dcopy)(&rowsOut, lpDblOut.front()->get(), &one, pDblOut->get() + pos, &one);
755 lpDblOut.pop_front();
757 out.push_back(pDblOut);
761 types::Double* pDblHdOut = new types::Double(rworksize + iworksize, 1);
762 C2F(dcopy)(&rworksize, rwork, &one, pDblHdOut->get(), &one);
764 for (int i = 0; i < iworksize; i++)
766 pDblHdOut->set(rworksize + i, (double)iwork[i]);
769 out.push_back(pDblHdOut);
773 if (pDblAtol == NULL || pDblAtol->isScalar())
778 if (pDblRtol == NULL || pDblRtol->isScalar())
789 DifferentialEquation::removeDifferentialEquationFunctions();
791 return types::Function::OK;