scilab/modules/differential_equations/sci_gateway/c/Ex-daskr.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include "machine.h"
   4 #include "core_math.h"
   5
   6 extern void C2F(dgesv)(int*, int*, double*, int*, int*, double*, int*, int*);
   7 extern double C2F(dlamch)(const char*);
   8 typedef void (*resfunc)(double*, double*, double*, double*, int*, double*, int*);
   9
  10 void C2F(pjac1)(resfunc res, int *ires, int *nequations, double *tOld, double *actual, double *actualP,
  11                 double *rewt, double *savr, double *wk, double *h, double *cj, double *wp, int *iwp,
  12                 int *ier, double *rpar, int *ipar)
  13 {
  14     int i = 0;
  15     int j = 0;
  16     int nrow = 0;
  17     double tx = 0;
  18     double del = 0;
  19     double delinv = 0;
  20     double ysave = 0;
  21     double ypsave = 0;
  22     double* e = NULL;
  23
  24     int neq = *nequations;
  25     double SQuround = sqrt(C2F(dlamch)("P"));
  26
  27     tx = *tOld;
  28
  29     e = (double*) calloc(neq, sizeof(double));
  30
  31     for (i = 0 ; i < neq ; ++i)
  32     {
  33         del = Max(SQuround * Max(fabs(actual[i]), fabs(*h * actualP[i])), 1. / rewt[i]);
  34         del *= (*h * actualP[i] >= 0) ? 1 : -1;
  35         ysave   =  actual[i];
  36         ypsave  =  actualP[i];
  37         actual[i]  += del;
  38         actualP[i] += *cj * del;
  39         res(&tx, actual, actualP, e, ires, rpar, ipar);
  40
  41         if (*ires < 0)
  42         {
  43             *ier = -1;
  44             free(e);
  45             return;
  46         }
  47
  48         delinv = 1. / del;
  49         for (j = 0 ; j < neq ; j++)
  50         {
  51             wp[nrow + j] = (e[j] - savr[j]) * delinv;
  52             /* NaN test */
  53             if (ISNAN(wp[nrow + j]))
  54             {
  55                 *ier = -1;
  56                 free (e);
  57                 return;
  58             }
  59             iwp[nrow + j] = i + 1;
  60             iwp[nrow + j + neq * neq] = j + 1;
  61         }
  62         nrow       += neq;
  63         actual[i]  =  ysave;
  64         actualP[i] =  ypsave;
  65     }
  66
  67     free(e);
  68 }
  69
  70 void C2F(psol1)(int *nequations, double *tOld, double *actual, double *actualP,
  71                 double *savr, double *wk, double *cj, double *wght, double *wp,
  72                 int *iwp, double *b, double *eplin, int *ier, double *dummy1, int *dummy2)
  73 {
  74     int i = 0;
  75     int nColB = 1;
  76     int info = 0;
  77     int *ipiv = NULL;
  78
  79     ipiv = (int *) malloc(*nequations * sizeof(int));
  80
  81     C2F(dgesv) (nequations, &nColB, wp, nequations, ipiv, b, nequations, &info);
  82
  83     if (info == 0)
  84     {
  85         /* NaN test */
  86         for (i = 0; i < *nequations; ++i)
  87         {
  88             if (ISNAN(b[i]))
  89             {
  90                 *ier = -1;
  91                 free (ipiv);
  92                 return;
  93             }
  94         }
  95     }
  96     else
  97     {
  98         *ier = -1;
  99     }
 100
 101     free (ipiv);
 102 }