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 int C2F(dgesv)(int*, int*, double*, int*, int*, double*, int*, int*);
   7 extern double C2F(dlamch)(char*);
   8 typedef void (*resfunc)(double*, double*, double*, double*, int*, double*, int*);
   9
  10 int 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     int neq = *nequations;
  24     double SQuround = sqrt(C2F(dlamch)("P"));
  25
  26     tx = *tOld;
  27
  28     e = (double*)calloc(neq, sizeof(double));
  29
  30     for (i = 0 ; i < neq ; ++i)
  31     {
  32         del = Max(SQuround * Max(fabs(actual[i]), fabs(*h * actualP[i])), 1. / rewt[i]);
  33         del *= (*h * actualP[i] >= 0) ? 1 : -1;
  34         ysave   =  actual[i];
  35         ypsave  =  actualP[i];
  36         actual[i]  += del;
  37         actualP[i] += *cj * del;
  38         res(&tx, actual, actualP, e, ires, rpar, ipar);
  39
  40         if (ires < 0)
  41         {
  42             return 0;
  43         }
  44
  45         delinv = 1. / del;
  46         for (j = 0 ; j < neq ; j++)
  47         {
  48             wp[nrow + j] = (e[j] - savr[j]) * delinv;
  49             iwp[nrow + j] = i + 1;
  50             iwp[nrow + j + neq * neq] = j + 1;
  51         }
  52         nrow       += neq;
  53         actual[i]  =  ysave;
  54         actualP[i] =  ypsave;
  55     }
  56
  57     free(e);
  58
  59     return (*ier);
  60 }
  61
  62 int C2F(psol1)(int *nequations, double *tOld, double *actual, double *actualP,
  63                double *savr, double *wk, double *cj, double *wght, double *wp,
  64                int *iwp, double *b, double *eplin, int *ier, double *dummy1, int *dummy2)
  65 {
  66     int nColB = 1;
  67     int info = 0;
  68     int *ipiv = NULL;
  69
  70     ipiv = (int *) malloc(*nequations * sizeof(int));
  71
  72     C2F(dgesv) (nequations, &nColB, wp, nequations, ipiv, b, nequations, &info);
  73
  74     free(ipiv);
  75     return 0;
  76 }
  77