Differential_equations: daskr subroutines error control 47/11947/5
Paul BIGNIER [Fri, 5 Jul 2013 09:49:33 +0000 (11:49 +0200)]
Improved the error control in daskr's example routines.

Reminder: these functions are used as demos, but the user can pass its own functions

Change-Id: I9e4794a1b38596630fb1d70f6439f9159e31abc0

scilab/modules/differential_equations/sci_gateway/c/Ex-daskr.c
scilab/modules/differential_equations/tests/unit_tests/daskr.dia.ref
scilab/modules/differential_equations/tests/unit_tests/daskr.tst

index 01ae0fd..11f9f7c 100644 (file)
@@ -3,13 +3,13 @@
 #include "machine.h"
 #include "core_math.h"
 
-extern int C2F(dgesv)(int*, int*, double*, int*, int*, double*, int*, int*);
-extern double C2F(dlamch)(char*);
+extern void C2F(dgesv)(int*, int*, double*, int*, int*, double*, int*, int*);
+extern double C2F(dlamch)(const char*);
 typedef void (*resfunc)(double*, double*, double*, double*, int*, double*, int*);
 
-int C2F(pjac1)(resfunc res, int *ires, int *nequations, double *tOld, double *actual, double *actualP,
-               double *rewt, double *savr, double *wk, double *h, double *cj, double *wp, int *iwp,
-               int *ier, double *rpar, int *ipar)
+void C2F(pjac1)(resfunc res, int *ires, int *nequations, double *tOld, double *actual, double *actualP,
+                double *rewt, double *savr, double *wk, double *h, double *cj, double *wp, int *iwp,
+                int *ier, double *rpar, int *ipar)
 {
     int i = 0;
     int j = 0;
@@ -20,12 +20,13 @@ int C2F(pjac1)(resfunc res, int *ires, int *nequations, double *tOld, double *ac
     double ysave = 0;
     double ypsave = 0;
     double* e = NULL;
+
     int neq = *nequations;
     double SQuround = sqrt(C2F(dlamch)("P"));
 
     tx = *tOld;
 
-    e = (double*)calloc(neq, sizeof(double));
+    e = (double*) calloc(neq, sizeof(double));
 
     for (i = 0 ; i < neq ; ++i)
     {
@@ -37,15 +38,24 @@ int C2F(pjac1)(resfunc res, int *ires, int *nequations, double *tOld, double *ac
         actualP[i] += *cj * del;
         res(&tx, actual, actualP, e, ires, rpar, ipar);
 
-        if (ires < 0)
+        if (*ires < 0)
         {
-            return 0;
+            *ier = -1;
+            free(e);
+            return;
         }
 
         delinv = 1. / del;
         for (j = 0 ; j < neq ; j++)
         {
             wp[nrow + j] = (e[j] - savr[j]) * delinv;
+            /* NaN test */
+            if (ISNAN(wp[nrow + j]))
+            {
+                *ier = -1;
+                free (e);
+                return;
+            }
             iwp[nrow + j] = i + 1;
             iwp[nrow + j + neq * neq] = j + 1;
         }
@@ -55,14 +65,13 @@ int C2F(pjac1)(resfunc res, int *ires, int *nequations, double *tOld, double *ac
     }
 
     free(e);
-
-    return (*ier);
 }
 
-int C2F(psol1)(int *nequations, double *tOld, double *actual, double *actualP,
-               double *savr, double *wk, double *cj, double *wght, double *wp,
-               int *iwp, double *b, double *eplin, int *ier, double *dummy1, int *dummy2)
+void C2F(psol1)(int *nequations, double *tOld, double *actual, double *actualP,
+                double *savr, double *wk, double *cj, double *wght, double *wp,
+                int *iwp, double *b, double *eplin, int *ier, double *dummy1, int *dummy2)
 {
+    int i = 0;
     int nColB = 1;
     int info = 0;
     int *ipiv = NULL;
@@ -71,7 +80,23 @@ int C2F(psol1)(int *nequations, double *tOld, double *actual, double *actualP,
 
     C2F(dgesv) (nequations, &nColB, wp, nequations, ipiv, b, nequations, &info);
 
-    free(ipiv);
-    return 0;
-}
+    if (info == 0)
+    {
+        /* NaN test */
+        for (i = 0; i < *nequations; ++i)
+        {
+            if (ISNAN(b[i]))
+            {
+                *ier = -1;
+                free (ipiv);
+                return;
+            }
+        }
+    }
+    else
+    {
+        *ier = -1;
+    }
 
+    free (ipiv);
+}
index edee5f2..10fef2a 100644 (file)
@@ -8,7 +8,6 @@
 // The terms are also available at
 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 // =============================================================================
-// Run with test_run('differential_equations', 'daskr', ['no_check_error_output']);
 //C-----------------------------------------------------------------------
 //C First problem.
 //C The initial value problem is..
@@ -45,15 +44,6 @@ y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
 if abs(nn(1)-2.53)>0.001 then bugmes();quit;end
 // Same problem, but using macros for the preconditioner evaluation and application functions 'pjac' and 'psol'
 // pjac uses the macro res1 defined above.
-function [r, ier] = psol(wp, iwp, b)
-    ier = 0;
-    //Compute the LU factorization of R.
-    sp = sparse(iwp, wp);
-    [h, rk] = lufact(sp);
-    //Solve the system LU*X = b
-    r = lusolve(h, b);
-    ludel(h);
-endfunction
 function [wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr)
     ires = 0;
     SQuround = 1.490D-08;
@@ -70,18 +60,34 @@ function [wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr)
         y(i) = y(i) + del;
         ydot(i) = ydot(i) + cj*del;
         [e ires]=res1(tx, y, ydot);
-        if ires < 0 then return; end
+        if ires < 0 then
+            ires = -1;
+            return;
+        end
         delinv = 1/del;
         for j=1:neq
             wp(nrow+j) = delinv*(e(j)-savr(j));
-            iwp(nrow+j,1) = i;
-            iwp(nrow+j,2) = j;
+            if isnan(wp(nrow+j)) then
+                ires = -1;
+                return;
+            end
+            iwp(nrow+j, 1) = i;
+            iwp(nrow+j, 2) = j;
         end
         nrow = nrow + neq;
         y(i) = ysave;
         ydot(i) = ypsave;
     end
 endfunction
+function [r, ier] = psol(wp, iwp, b)
+    ier = 0;
+    //Compute the LU factorization of R.
+    sp = sparse(iwp, wp);
+    [h, rk] = lufact(sp);
+    //Solve the system LU*X = b
+    r = lusolve(h, b);
+    ludel(h);
+endfunction
 y0=1;t=2:6;t0=1;y0d=3;
 info=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
 atol=1.d-6;rtol=0;ng=2;
index 2e91c81..c11ea98 100644 (file)
@@ -9,8 +9,6 @@
 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 // =============================================================================
 
-// Run with test_run('differential_equations', 'daskr', ['no_check_error_output']);
-
 //C-----------------------------------------------------------------------
 //C First problem.
 //C The initial value problem is..
@@ -50,15 +48,6 @@ if abs(nn(1)-2.53)>0.001 then pause,end
 
 // Same problem, but using macros for the preconditioner evaluation and application functions 'pjac' and 'psol'
 // pjac uses the macro res1 defined above.
-function [r, ier] = psol(wp, iwp, b)
-    ier = 0;
-    //Compute the LU factorization of R.
-    sp = sparse(iwp, wp);
-    [h, rk] = lufact(sp);
-    //Solve the system LU*X = b
-    r = lusolve(h, b);
-    ludel(h);
-endfunction
 function [wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr)
     ires = 0;
     SQuround = 1.490D-08;
@@ -75,18 +64,34 @@ function [wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr)
         y(i) = y(i) + del;
         ydot(i) = ydot(i) + cj*del;
         [e ires]=res1(tx, y, ydot);
-        if ires < 0 then return; end
+        if ires < 0 then
+            ires = -1;
+            return;
+        end
         delinv = 1/del;
         for j=1:neq
             wp(nrow+j) = delinv*(e(j)-savr(j));
-            iwp(nrow+j,1) = i;
-            iwp(nrow+j,2) = j;
+            if isnan(wp(nrow+j)) then
+                ires = -1;
+                return;
+            end
+            iwp(nrow+j, 1) = i;
+            iwp(nrow+j, 2) = j;
         end
         nrow = nrow + neq;
         y(i) = ysave;
         ydot(i) = ypsave;
     end
 endfunction
+function [r, ier] = psol(wp, iwp, b)
+    ier = 0;
+    //Compute the LU factorization of R.
+    sp = sparse(iwp, wp);
+    [h, rk] = lufact(sp);
+    //Solve the system LU*X = b
+    r = lusolve(h, b);
+    ludel(h);
+endfunction
 
 y0=1;t=2:6;t0=1;y0d=3;
 info=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);