4 #include "mkl_lapacke.h"
13 ImplicitMKL(
int my_neq,
int my_nt,
double my_ti,
double my_tf,
double my_dt) {
21 void rhs(
double t,
double *u,
double *f) {
28 for (
int i =0; i<neq; i++) {
33 void step(
double t,
double * u,
double * unew) {
43 double NEWTON_TOL = 1.0e-14;
44 int NEWTON_MAXSTEP = 1000;
50 stepsize =
new double[neq];
51 uguess =
new double[neq];
53 J =
new double[neq*neq];
56 for (
int i=0; i<neq; i++) {
63 int * pivot =
new int[neq];
68 newt(tnew,u,uguess,stepsize);
74 LAPACKE_dgesv(LAPACK_ROW_MAJOR, (lapack_int) neq,
75 (lapack_int) 1, J, neq,pivot,stepsize,1);
79 for (
int i = 0; i<neq; i++) {
80 uguess[i] = uguess[i] - stepsize[i];
81 if (std::abs(stepsize[i])>maxstepsize) {
82 maxstepsize = std::abs(stepsize[i]);
87 if ( maxstepsize < NEWTON_TOL) {
93 if (counter > NEWTON_MAXSTEP) {
94 fprintf(stderr,
"max Newton iterations reached\n");
99 for (
int i=0; i<neq; i++) {
111 void newt(
double t,
double *uprev,
double *uguess,
123 for (
int i =0; i<neq; i++) {
124 g[i] = uguess[i]-uprev[i]-dt*g[i];
128 void jac(
double t,
double *u,
double *J){
145 u1 =
new double[neq];
146 f1 =
new double[neq];
152 for (
int i = 0; i<neq; i++) {
153 for (
int j = 0; j<neq; j++) {
160 for (
int j = 0; j<neq; j++) {
161 J[j*neq+i] = -dt*(f1[j]-f[j])/d;
163 J[i*neq+i] = 1.0+J[i*neq+i];
179 #endif // _IMPLICIT_H_
ImplicitMKL(int my_neq, int my_nt, double my_ti, double my_tf, double my_dt)
void newt(double t, double *uprev, double *uguess, double *g)
header file containing explanation of functions for the RIDC integrator
void rhs(double t, double *u, double *f)
void newt(double t, double *uprev, double *Kguess, double *g, PARAMETER param, BUTCHER rk)
void jac(double t, double *uprev, double *Kguess, double *J, PARAMETER param, BUTCHER rk)
void step(double t, double *u, double *unew)
void rhs(double t, double *u, PARAMETER param, double *f)
void jac(double t, double *u, double *J)