Bug #12418 fixed again: using bvode() with continuation led to an error 38/21538/4
Clément DAVID [Fri, 17 Jul 2020 13:12:06 +0000 (15:12 +0200)]
Change-Id: I0f3d19064db99e052a8273bfc233bee3b439f2de

scilab/CHANGES.md
scilab/modules/differential_equations/sci_gateway/cpp/sci_bvode.cpp
scilab/modules/differential_equations/tests/nonreg_tests/bug_12418.dia.ref
scilab/modules/differential_equations/tests/nonreg_tests/bug_12418.tst

index f4bf824..eb04414 100644 (file)
@@ -279,6 +279,7 @@ Bug Fixes
 * [#8378](https://bugzilla.scilab.org/8378): Datatip `ContextMenu => Delete last datatip` was useless.
 * [#9909](https://bugzilla.scilab.org/9909): In the help browser, add a way to open the online version of the current page.
 * [#10476](https://bugzilla.scilab.org/10476): From `browsevar`, displaying the content of lists, structures, cells, or other custom tlists or mlists was not possible.
+* [#12418](https://bugzilla.scilab.org/12418): Using bvode() with "continuation", i.e. `ipar(9) > 1` led to an error.
 * [#12516](https://bugzilla.scilab.org/12516): From `browsevar`, clicking on any graphical handle did not edit its figure with `ged`.
 * [#12532](https://bugzilla.scilab.org/12532): From `browsevar`, clicking on any function did not edit it with `edit`. The content of libraries could not be displayed either.
 * [#12719](https://bugzilla.scilab.org/12719): `A(%s)` gave the same result as `A($)`.
index 794d379..ab8c5ad 100644 (file)
@@ -22,6 +22,9 @@
 #include "callable.hxx"
 #include "differentialequationfunctions.hxx"
 
+#include <memory>
+#include <vector>
+
 extern "C"
 {
 #include "sci_malloc.h"
@@ -30,6 +33,13 @@ extern "C"
 #include "scifunctions.h"
 #include "sciprint.h"
 }
+/*--------------------------------------------------------------------------*/
+
+namespace {
+    // rwork / fspace is declared as a Fortran COMMON to call bvode with the same iguess/ipar(9) space 
+    static std::shared_ptr< std::vector<double> > rwork;
+    static std::shared_ptr< std::vector<int> > iwork;
+}
 
 /*--------------------------------------------------------------------------*/
 types::Function::ReturnValue sci_bvode(types::typed_list &in, int _iRetCount, types::typed_list &out)
@@ -547,8 +557,12 @@ types::Function::ReturnValue sci_bvode(types::typed_list &in, int _iRetCount, ty
     }
 
     // *** Perform operation. ***
-    double* rwork   = (double*)MALLOC(ipar[4] * sizeof(double));
-    int* iwork      = (int*)MALLOC(ipar[5] * sizeof(int));
+
+    // COMMON rwork and iwork (fspace and ispace)
+    if (!rwork || rwork->size() != ipar[4])
+        rwork = std::make_shared< std::vector<double> >(ipar[4]);
+    if (!iwork || iwork->size() != ipar[5])
+        iwork = std::make_shared< std::vector<int> >(ipar[5]);
     int* ltol       = (int*)MALLOC(pDblLtol->getSize() * sizeof(int));
 
     for (int i = 0; i < pDblLtol->getSize(); i++)
@@ -556,9 +570,15 @@ types::Function::ReturnValue sci_bvode(types::typed_list &in, int _iRetCount, ty
         ltol[i] = (int)pDblLtol->get(i);
     }
 
+    // for ipar(9) = 2 or 3 set ipar(3) = ispace(1) ( = the size of the previous mesh ).
+    if (ipar[8] == 2 || ipar[8] == 3)
+    {
+        ipar[2] = iwork->data()[0];
+    }
+
     try
     {
-        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);
+        C2F(colnew)(&ncomp, M, &aleft, &aright, pDblZeta->get(), ipar, ltol, pDblTol->get(), pDblFixpnt->get(), iwork->data(), rwork->data(), &iflag, bvode_fsub, bvode_dfsub, bvode_gsub, bvode_dgsub, bvode_guess);
     }
     catch (ast::InternalError &ie)
     {
@@ -568,8 +588,6 @@ types::Function::ReturnValue sci_bvode(types::typed_list &in, int _iRetCount, ty
 
     if (bCatch)
     {
-        FREE(iwork);
-        FREE(rwork);
         FREE(M);
         FREE(ltol);
         DifferentialEquation::removeDifferentialEquationFunctions();
@@ -599,8 +617,6 @@ types::Function::ReturnValue sci_bvode(types::typed_list &in, int _iRetCount, ty
             Scierror(999, _("%s: There is an input data error.\n"), "bvode");
         }
 
-        FREE(iwork);
-        FREE(rwork);
         FREE(M);
         FREE(ltol);
         DifferentialEquation::removeDifferentialEquationFunctions();
@@ -613,13 +629,11 @@ types::Function::ReturnValue sci_bvode(types::typed_list &in, int _iRetCount, ty
     for (int i = 0; i < pDblXpts->getSize(); i++)
     {
         double val = pDblXpts->get(i);
-        C2F(appsln)(&val, &res[i * sumM], rwork, iwork);
+        C2F(appsln)(&val, &res[i * sumM], rwork->data(), iwork->data());
     }
 
     pDblRes->set(res);
 
-    FREE(iwork);
-    FREE(rwork);
     FREE(M);
     FREE(ltol);
     FREE(res);
index 2750b48..fe0019c 100644 (file)
@@ -43,31 +43,31 @@ ipar(4) = 4;
 ltol = [1, 2, 3, 4];
 tol = [1d-5, 1d-5, 1d-5, 1d-5];
 function [f]=fsub(x,z)
-       f = [ -.5*(1/x+1/(x-1)+1/(x-y))*z(2) +...
-     z(1) * ((v-z(3)-z(4))/x + z(3)/(x-1) + z(4)/(x-y)),...
-       0,0];
+    f = [ -.5*(1/x+1/(x-1)+1/(x-y))*z(2) +...
+    z(1) * ((v-z(3)-z(4))/x + z(3)/(x-1) + z(4)/(x-y)),...
+    0,0];
 endfunction
 function [df] = dfsub(x,z)
-       df = [(v-z(3)-z(4))/x + z(3)/(x-1) + z(4)/(x-y),...
-       -.5*(1/x+1/(x-1)+1/(x-y)),z(1)/(x*(x-1)),z(1)*y/(x*(x-y));...
-       0,0,0,0;0,0,0,0];
+    df = [(v-z(3)-z(4))/x + z(3)/(x-1) + z(4)/(x-y),...
+    -.5*(1/x+1/(x-1)+1/(x-y)),z(1)/(x*(x-1)),z(1)*y/(x*(x-y));...
+    0,0,0,0;0,0,0,0];
 endfunction
 // Boundary conditions
 function [g]=gsub(i,z)
-       select i
-       case 1, g = z(2) - 2*z(1)*(v-z(3)-z(4))
-       case 2, g = z(2) - 2*z(1)*z(3)
-       case 3, g = z(1)-1.
-       case 4, g = z(2) - 2*z(1)*z(4)
-       end
+    select i
+    case 1, g = z(2) - 2*z(1)*(v-z(3)-z(4))
+    case 2, g = z(2) - 2*z(1)*z(3)
+    case 3, g = z(1)-1.
+    case 4, g = z(2) - 2*z(1)*z(4)
+    end
 endfunction
 function [dg]=dgsub(i,z)
-       select i
-       case 1, dg = [-2*(v-z(3)-z(4)),1.,2*z(1),2*z(1)]
-       case 2, dg = [-2*z(3),1.,-2*z(1),0]
-       case 3, dg = [1,0,0,0]
-       case 4, dg = [-2*z(4),1.,0,-2*z(1)]
-       end
+    select i
+    case 1, dg = [-2*(v-z(3)-z(4)),1.,2*z(1),2*z(1)]
+    case 2, dg = [-2*z(3),1.,-2*z(1),0]
+    case 3, dg = [1,0,0,0]
+    case 4, dg = [-2*z(4),1.,0,-2*z(1)]
+    end
 endfunction
 // Start computation
 // Locations of side conditions, sorted
@@ -78,17 +78,17 @@ aright = y;
 valv = [linspace(0,.9,10) logspace(0,2,30)];
 res = [linspace(0,.99,100) linspace(1,y,101)];
 function [z,dmval]=guess(x)
-        z=[2*x-1, 2., 1., 1/(2*y-1)]
-        dmval=[0,0,0]
+    z=[2*x-1, 2., 1., 1/(2*y-1)]
+    dmval=[0,0,0]
 endfunction
 // First execution has ipar(9)=1 and uses the guess
 // Subsequent executions have ipar(9)=3 and use continuation. This is
 // run in tight closed loop to not disturb the stack
 for i=1:40
-v=valv(i);
-sol=bvode(res,ncomp,m,aleft,aright,zeta,ipar,ltol,tol,fixpnt,...
- fsub,dfsub,gsub,dgsub,guess);
-eigens(:,i)=[v;sol(3,101);sol(4,101)];  // c_2 and c_3 are constant!
-ipar(9)=3;
+    v=valv(i);
+    sol=bvode(res,ncomp,m,aleft,aright,zeta,ipar,ltol,tol,fixpnt,...
+    fsub,dfsub,gsub,dgsub,guess);
+    eigens(:,i)=[v;sol(3,101);sol(4,101)];  // c_2 and c_3 are constant!
+    ipar(9)=3;
 end
 assert_checkalmostequal(eigens, [0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,1.1721023,1.3738238,1.610262,1.8873918,2.2122163,2.5929438,3.0391954,3.5622479,4.1753189,4.8939009,5.7361525,6.7233575,7.8804628,9.2367086,10.826367,12.68961,14.873521,17.433288,20.433597,23.950266,28.072162,32.903446,38.566204,45.203537,52.983169,62.101694,72.789538,85.316785,100;0.8590043,0.8923732,0.9270092,0.9629650,1.0002848,1.039002,1.0791385,1.1207035,1.1636934,1.2080921,1.2538713,1.3357726,1.4364687,1.5602838,1.712184,1.8976546,2.1225064,2.392618,2.713585,3.0901836,3.5255313,4.0199156,4.569725,5.1676944,5.8058304,6.4799746,7.1921595,7.9490495,8.7589833,9.6301157,10.56993,11.585457,12.683653,13.871706,15.157239,16.548439,18.054127,19.683839,21.447881,23.357401;0.3373662,0.3544495,0.3714346,0.3883257,0.4051247,0.4218310,0.4384421,0.4549534,0.4713585,0.4876500,0.5038194,0.5313357,0.5630336,0.5993442,0.6406497,0.6872575,0.7393884,0.7971803,0.8606997,0.9299459,1.0048381,1.0851888,1.1707035,1.2610735,1.3561975,1.4564017,1.5624488,1.6753378,1.7961064,1.9257493,2.0652244, 2.2154891,2.3775328,2.5524022,2.7412151,2.9451754,3.1655779,3.4038197,3.6614094,3.9399763], 10^-6 );
index 52b5049..69a721b 100644 (file)
@@ -14,8 +14,7 @@
 //
 // <-- Short Description -->
 // Continuation is incorrectly supported in bvode
-//
-// <-- NOT FIXED -->
+
 
 
 y= 1.9d0;