2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2011 - DIGITEO - Cedric DELAMARRE
5 * Copyright (C) 2012 - 2016 - Scilab Enterprises
7 * This file is hereby licensed under the terms of the GNU GPL v2.0,
8 * pursuant to article 5.3.4 of the CeCILL v.2.1.
9 * This file was originally licensed under the terms of the CeCILL v2.1,
10 * and continues to be available under such terms.
11 * For more information, see the COPYING file which you should have received
12 * along with this program.
15 /*--------------------------------------------------------------------------*/
17 #include "differential_equations_gw.hxx"
18 #include "function.hxx"
22 #include "callable.hxx"
23 #include "differentialequationfunctions.hxx"
27 #include "sci_malloc.h"
28 #include "localization.h"
30 #include "scifunctions.h"
34 /*--------------------------------------------------------------------------*/
35 types::Function::ReturnValue sci_bvode(types::typed_list &in, int _iRetCount, types::typed_list &out)
40 types::Double* pDblXpts = NULL;
41 types::Double* pDblM = NULL;
42 types::Double* pDblZeta = NULL;
43 types::Double* pDblIpar = NULL;
44 types::Double* pDblLtol = NULL;
45 types::Double* pDblTol = NULL;
46 types::Double* pDblFixpnt = NULL;
56 // error message catched
57 std::wostringstream os;
60 // *** check the minimal number of input args. ***
63 Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "bvode", 15);
64 return types::Function::Error;
67 // *** check number of output args ***
70 Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "bvode", 1);
71 return types::Function::Error;
74 // *** check type of input args and get it. ***
77 if (in[iPos]->isDouble() == false)
79 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "bvode", iPos + 1);
80 return types::Function::Error;
83 pDblXpts = in[iPos]->getAs<types::Double>();
87 if (in[iPos]->isDouble() == false)
89 Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "bvode", iPos + 1);
90 return types::Function::Error;
93 types::Double* pDblN = in[iPos]->getAs<types::Double>();
95 if (pDblN->isScalar() == false)
97 Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "bvode", iPos + 1);
98 return types::Function::Error;
101 ncomp = (int)pDblN->get(0);
105 Scierror(999, _("%s: Wrong value for input argument #%d: Value at most 20 expected.\n"), "bvode", iPos + 1);
106 return types::Function::Error;
111 if (in[iPos]->isDouble() == false)
113 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "bvode", iPos + 1);
114 return types::Function::Error;
117 pDblM = in[iPos]->getAs<types::Double>();
119 if (pDblM->getSize() != ncomp)
121 Scierror(999, _("%s: Wrong size for input argument #%d: A vector of size %d (N) expected.\n"), "bvode", iPos + 1, ncomp);
122 return types::Function::Error;
125 int* M = (int*)MALLOC(pDblM->getSize() * sizeof(int));
126 for (int i = 0; i < pDblM->getSize(); i++)
128 M[i] = (int)pDblM->get(i);
129 sumM += (int)pDblM->get(i);
130 maxM = std::max(maxM, (int)pDblM->get(i));
135 Scierror(999, _("%s: Wrong value for input argument #%d: Sum of m (%d) must be less than 40.\n"), "bvode", iPos + 1, sumM);
137 return types::Function::Error;
143 if (in[iPos]->isDouble() == false)
145 Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "bvode", iPos + 1);
147 return types::Function::Error;
150 types::Double* pDblXLow = in[iPos]->getAs<types::Double>();
152 if (pDblXLow->isScalar() == false)
154 Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "bvode", iPos + 1);
156 return types::Function::Error;
159 aleft = pDblXLow->get(0);
163 if (in[iPos]->isDouble() == false)
165 Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "bvode", iPos + 1);
167 return types::Function::Error;
170 types::Double* pDblXUp = in[iPos]->getAs<types::Double>();
172 if (pDblXUp->isScalar() == false)
174 Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "bvode", iPos + 1);
176 return types::Function::Error;
179 aright = pDblXUp->get(0);
183 if (in[iPos]->isDouble() == false)
185 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "bvode", iPos + 1);
187 return types::Function::Error;
190 pDblZeta = in[iPos]->getAs<types::Double>();
192 for (int i = 0; i < pDblZeta->getSize() - 1; i++)
194 if (pDblZeta->get(i) > pDblZeta->get(i + 1))
196 Scierror(999, _("%s: Wrong value for input argument #%d: zeta(j) lower or equal to zeta(j+1) expected.\n"), "bvode", iPos + 1);
198 return types::Function::Error;
205 if (in[iPos]->isDouble() == false)
207 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "bvode", iPos + 1);
209 return types::Function::Error;
212 pDblIpar = in[iPos]->getAs<types::Double>();
214 if (pDblIpar->getSize() != 11)
216 Scierror(999, _("%s: Wrong size for input argument #%d: A vector of size 11 expected.\n"), "bvode", iPos + 1);
218 return types::Function::Error;
221 for (int i = 0; i < 11; i++)
223 ipar[i] = (int)pDblIpar->get(i);
226 // check ipar arguments
228 // ipar(1) : 0 => linear 1 => nonlinear
229 if (ipar[0] != 0 && ipar[0] != 1)
231 Scierror(999, _("%s: Wrong value for input argument #%d: Wrong value for element %s: %d or %d expected.\n"), "bvode", iPos + 1, 1, 0, 1);
233 return types::Function::Error;
236 // ipar(2)=0 then collpnt is set to max(max(m(j))+1, 5-max(m(j))).
239 ipar[1] = std::max(maxM + 1, 5 - maxM);
242 if (ipar[1] < 0 || ipar[1] > 7)
244 Scierror(999, _("%s: Wrong value for input argument #%d: Wrong value for element %s: %d to %d expected.\n"), "bvode", iPos + 1, 2, 0, 7);
246 return types::Function::Error;
256 if (ipar[3] < 0 || ipar[3] > sumM)
258 Scierror(999, _("%s: Wrong value for input argument #%d: Wrong value for element %s: %d to sum of m (=%d) expected.\n"), "bvode", iPos + 1, 4, 0, sumM);
260 return types::Function::Error;
263 // bvode: ipar(5)>=nmax*nsizef where
264 // nsizef=4 + 3*M + (5+collpnt*N)*(collpnt*N+M) + (2*M-nrec)*2*M
266 int kdm = ipar[1] * ncomp;
268 for (int i = 0; i < sumM; i++)
270 int ib = sumM + 1 - i;
271 if (pDblZeta->getSize() >= ib && pDblZeta->get(ib - 1) >= aright)
277 int nsizef = 4 + 3 * sumM + (5 + kdm) * (kdm + sumM) + (2 * sumM - nrec) * 2 * sumM;
278 if (ipar[4] < nmax * nsizef)
280 Scierror(999, _("%s: Wrong value for input argument #%d: Wrong value for element %d.\n"), "bvode", iPos + 1, 5);
282 return types::Function::Error;
285 // ipar(6)>=nmax*nsizei where nsizei= 3 + collpnt*N + M.
286 int nsizei = 3 + kdm + sumM;
287 if (ipar[5] < nmax * nsizei)
289 Scierror(999, _("%s: Wrong value for input argument #%d: Wrong value for element %d.\n"), "bvode", iPos + 1, 6);
291 return types::Function::Error;
294 // ipar(7) output control
295 if (ipar[6] < -1 || ipar[6] > 1)
297 Scierror(999, _("%s: Wrong value for input argument #%d: Wrong value for element %d: %d to %d expected.\n"), "bvode", iPos + 1, 7, -1, 1);
299 return types::Function::Error;
303 if (ipar[7] < 0 || ipar[7] > 2)
305 Scierror(999, _("%s: Wrong value for input argument #%d: Wrong value for element %d: %d to %d expected.\n"), "bvode", iPos + 1, 8, 0, 2);
307 return types::Function::Error;
311 if (ipar[8] < 0 || ipar[8] > 4)
313 Scierror(999, _("%s: Wrong value for input argument #%d: Wrong value for element %d: %d to %d expected.\n"), "bvode", iPos + 1, 9, 0, 4);
315 return types::Function::Error;
319 if (ipar[9] < 0 || ipar[9] > 2)
321 Scierror(999, _("%s: Wrong value for input argument #%d: Wrong value for element %d: %d to %d expected.\n"), "bvode", iPos + 1, 9, 0, 2);
323 return types::Function::Error;
330 if (in[iPos]->isDouble() == false)
332 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "bvode", iPos + 1);
334 return types::Function::Error;
337 pDblLtol = in[iPos]->getAs<types::Double>();
339 if (pDblLtol->getSize() != ipar[3])
341 Scierror(999, _("%s: Wrong size for input argument #%d: An array of size %d (ipar(4)) expected.\n"), "bvode", iPos + 1, ipar[3]);
343 return types::Function::Error;
346 // verify following cases :
347 // 1 <= ltol(1) < ltol(2) < ... < ltol(NTOL) <= M where M=sum(m)
348 // M is mstar and NTOL is the size of ltol
351 double* pdblLtol = pDblLtol->get();
352 if (pdblLtol[0] < 1 || pdblLtol[0] > sumM)
354 Scierror(999, _("%s: Wrong value for input argument #%d: %d to sum of m (=%d) expected.\n"), "bvode", iPos + 1, 1, sumM);
356 return types::Function::Error;
359 // ltol(1) < ltol(2) < ... < ltol(NTOL) <= M
360 for (int i = 1; i < pDblLtol->getSize(); i++)
362 if (pdblLtol[i - 1] >= pdblLtol[i] || pdblLtol[i - 1] > sumM)
364 Scierror(999, _("%s: Wrong value for input argument #%d: Bad value for ltol(%d): ltol(1) < ltol(2) < ... < ltol(NTOL) <= M (sum of m) expected.\n"), "bvode", iPos + 1, i);
366 return types::Function::Error;
373 if (in[iPos]->isDouble() == false)
375 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "bvode", iPos + 1);
377 return types::Function::Error;
380 pDblTol = in[iPos]->getAs<types::Double>();
382 if (pDblTol->getSize() != ipar[3])
384 Scierror(999, _("%s: Wrong size for input argument #%d: An array of size %d (ipar(4)) expected.\n"), "bvode", iPos + 1, ipar[3]);
386 return types::Function::Error;
392 if (in[iPos]->isDouble() == false)
394 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "bvode", iPos + 1);
396 return types::Function::Error;
399 pDblFixpnt = in[iPos]->getAs<types::Double>();
401 if (pDblFixpnt->getSize() != ipar[10] && ipar[10] != 0)
403 Scierror(999, _("%s: Wrong size for input argument #%d: An array of size %d (ipar(11)) expected.\n"), "bvode", iPos + 1, ipar[10]);
405 return types::Function::Error;
410 // functions: fsub,dfsub,gsub,dgsub,guess
411 DifferentialEquationFunctions deFunctionsManager(L"bvode");
412 DifferentialEquation::addDifferentialEquationFunctions(&deFunctionsManager);
413 deFunctionsManager.setBvodeM(sumM);
414 deFunctionsManager.setBvodeN(ncomp);
416 for (int i = iPos; i < 15; i++)
418 if (in[i]->isCallable())
420 types::Callable* pCall = in[i]->getAs<types::Callable>();
423 deFunctionsManager.setFsubFunction(pCall);
425 else if (i == 11) // dfsub
427 deFunctionsManager.setDfsubFunction(pCall);
429 else if (i == 12) // gsub
431 deFunctionsManager.setGsubFunction(pCall);
433 else if (i == 13) // dgsub
435 deFunctionsManager.setDgsubFunction(pCall);
437 else if (i == 14 && ipar[8] == 1) // guess is needed only if ipar(9) == 1
439 deFunctionsManager.setGuessFunction(pCall);
442 else if (in[i]->isString())
445 types::String* pStr = in[i]->getAs<types::String>();
448 bOK = deFunctionsManager.setFsubFunction(pStr);
450 else if (i == 11) // dfsub
452 bOK = deFunctionsManager.setDfsubFunction(pStr);
454 else if (i == 12) // gsub
456 bOK = deFunctionsManager.setGsubFunction(pStr);
458 else if (i == 13) // dgsub
460 bOK = deFunctionsManager.setDgsubFunction(pStr);
462 else if (i == 14 && ipar[8] == 1) // guess is needed only if ipar(9) == 1
464 bOK = deFunctionsManager.setGuessFunction(pStr);
469 char* pst = wide_string_to_UTF8(pStr->get(0));
470 Scierror(50, _("%s: Subroutine not found: %s\n"), "bvode", pst);
472 DifferentialEquation::removeDifferentialEquationFunctions();
474 return types::Function::Error;
477 else if (in[i]->isList())
479 types::List* pList = in[i]->getAs<types::List>();
481 if (pList->getSize() == 0)
483 Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "bvode", i + 1, "(string empty)");
484 DifferentialEquation::removeDifferentialEquationFunctions();
486 return types::Function::Error;
489 if (pList->get(0)->isCallable())
493 deFunctionsManager.setFsubFunction(pList->get(0)->getAs<types::Callable>());
494 for (int iter = 1; iter < pList->getSize(); iter++)
496 deFunctionsManager.setFsubArgs(pList->get(iter)->getAs<types::InternalType>());
499 else if (i == 11) // dfsub
501 deFunctionsManager.setDfsubFunction(pList->get(0)->getAs<types::Callable>());
502 for (int iter = 1; iter < pList->getSize(); iter++)
504 deFunctionsManager.setDfsubArgs(pList->get(iter)->getAs<types::InternalType>());
507 else if (i == 12) // gsub
509 deFunctionsManager.setGsubFunction(pList->get(0)->getAs<types::Callable>());
510 for (int iter = 1; iter < pList->getSize(); iter++)
512 deFunctionsManager.setGsubArgs(pList->get(iter)->getAs<types::InternalType>());
515 else if (i == 13) // dgsub
517 deFunctionsManager.setDgsubFunction(pList->get(0)->getAs<types::Callable>());
518 for (int iter = 1; iter < pList->getSize(); iter++)
520 deFunctionsManager.setDgsubArgs(pList->get(iter)->getAs<types::InternalType>());
523 else if (i == 14 && ipar[8] == 1) // guess is needed only if ipar(9) == 1
525 deFunctionsManager.setGuessFunction(pList->get(0)->getAs<types::Callable>());
526 for (int iter = 1; iter < pList->getSize(); iter++)
528 deFunctionsManager.setGuessArgs(pList->get(iter)->getAs<types::InternalType>());
534 Scierror(999, _("%s: Wrong type for input argument #%d: The first argument in the list must be a Scilab function.\n"), "bvode", i + 1);
535 DifferentialEquation::removeDifferentialEquationFunctions();
537 return types::Function::Error;
542 Scierror(999, _("%s: Wrong type for input argument #%d: A function expected.\n"), "bvode", i + 1);
543 DifferentialEquation::removeDifferentialEquationFunctions();
545 return types::Function::Error;
549 // *** Perform operation. ***
550 double* rwork = (double*)MALLOC(ipar[4] * sizeof(double));
551 int* iwork = (int*)MALLOC(ipar[5] * sizeof(int));
552 int* ltol = (int*)MALLOC(pDblLtol->getSize() * sizeof(int));
554 for (int i = 0; i < pDblLtol->getSize(); i++)
556 ltol[i] = (int)pDblLtol->get(i);
561 C2F(colnew)(&ncomp, M, &aleft, &aright, pDblZeta->get(), ipar, ltol, pDblTol->get(), pDblFixpnt->get(), iwork, rwork, &iflag, bvode_fsub, bvode_dfsub, bvode_gsub, bvode_dgsub, bvode_guess);
563 catch (ast::InternalError &ie)
565 os << ie.GetErrorMessage();
575 DifferentialEquation::removeDifferentialEquationFunctions();
577 wchar_t szError[bsiz];
578 os_swprintf(szError, bsiz, _W("%s: An error occurred in '%s' subroutine.\n").c_str(), "bvode", "bvode");
580 throw ast::InternalError(os.str());
587 Scierror(999, _("%s: The collocation matrix is singular.\n"), "bvode");
589 else if (iflag == -1)
591 Scierror(999, _("%s: The expected no. of subintervals exceeds storage specifications.\n"), "bvode");
593 else if (iflag == -2)
595 Scierror(999, _("%s: The nonlinear iteration has not converged.\n"), "bvode");
597 else if (iflag == -3)
599 Scierror(999, _("%s: There is an input data error.\n"), "bvode");
606 DifferentialEquation::removeDifferentialEquationFunctions();
608 return types::Function::Error;
611 types::Double* pDblRes = new types::Double(sumM, pDblXpts->getSize());
612 double* res = (double*)MALLOC(pDblXpts->getSize() * sumM * sizeof(double));
613 for (int i = 0; i < pDblXpts->getSize(); i++)
615 double val = pDblXpts->get(i);
616 C2F(appsln)(&val, &res[i * sumM], rwork, iwork);
626 DifferentialEquation::removeDifferentialEquationFunctions();
628 // *** Return result in Scilab. ***
629 out.push_back(pDblRes);
631 return types::Function::OK;
633 /*--------------------------------------------------------------------------*/