From: Paul BIGNIER Date: Fri, 5 Jul 2013 09:49:33 +0000 (+0200) Subject: Differential_equations: daskr subroutines error control X-Git-Tag: 5.5.0-beta-1~483 X-Git-Url: http://gitweb.scilab.org/?p=scilab.git;a=commitdiff_plain;h=ce063fa752769f2adb3d00e49c66c4ac4f58a99a Differential_equations: daskr subroutines error control 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 --- diff --git a/scilab/modules/differential_equations/sci_gateway/c/Ex-daskr.c b/scilab/modules/differential_equations/sci_gateway/c/Ex-daskr.c index 01ae0fd..11f9f7c 100644 --- a/scilab/modules/differential_equations/sci_gateway/c/Ex-daskr.c +++ b/scilab/modules/differential_equations/sci_gateway/c/Ex-daskr.c @@ -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); +} diff --git a/scilab/modules/differential_equations/tests/unit_tests/daskr.dia.ref b/scilab/modules/differential_equations/tests/unit_tests/daskr.dia.ref index edee5f2..10fef2a 100644 --- a/scilab/modules/differential_equations/tests/unit_tests/daskr.dia.ref +++ b/scilab/modules/differential_equations/tests/unit_tests/daskr.dia.ref @@ -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; diff --git a/scilab/modules/differential_equations/tests/unit_tests/daskr.tst b/scilab/modules/differential_equations/tests/unit_tests/daskr.tst index 2e91c81..c11ea98 100644 --- a/scilab/modules/differential_equations/tests/unit_tests/daskr.tst +++ b/scilab/modules/differential_equations/tests/unit_tests/daskr.tst @@ -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);