Differential_equations: daskr subroutines error control
[scilab.git] / scilab / modules / differential_equations / sci_gateway / c / Ex-daskr.c
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);
+}