# HG changeset patch
# User John W. Eaton
# Date 1483403074 18000
# Mon Jan 02 19:24:34 2017 -0500
# Node ID 20017c4feae68cac685dcb23e23855d2bdb0d43c
# Parent 0cf815201920c29033f738ca0283872673f961fe
F77_INT fixes for Octave 4.4.
diff --git a/src/odepkg_auxiliary_functions.h b/src/odepkg_auxiliary_functions.h
--- a/src/odepkg_auxiliary_functions.h
+++ b/src/odepkg_auxiliary_functions.h
@@ -19,6 +19,13 @@ along with this program; If not, see
#include "odepkg_auxiliary_functions.h"
-typedef octave_idx_type (*odepkg_ddaskr_restype)
+typedef F77_RET_T (*odepkg_ddaskr_restype)
(const double& T, const double* Y, const double* YPRIME,
- const double& CJ, double* DELTA, octave_idx_type& IRES,
- const double* RPAR, const octave_idx_type* IPAR);
+ const double& CJ, double* DELTA, F77_INT& IRES,
+ const double* RPAR, const F77_INT* IPAR);
-typedef octave_idx_type (*odepkg_ddaskr_jactype)
+typedef F77_RET_T (*odepkg_ddaskr_jactype)
(const double& T, const double* Y, const double* YPRIME,
double* PD, const double& CJ, const double* RPAR,
- const octave_idx_type* IPAR);
+ const F77_INT* IPAR);
-typedef octave_idx_type (*odepkg_ddaskr_psoltype)
- (const octave_idx_type& NEQ, const double& T, const double* Y,
+typedef F77_RET_T (*odepkg_ddaskr_psoltype)
+ (const F77_INT& NEQ, const double& T, const double* Y,
const double* YPRIME, const double* SAVR, const double* WK,
- const octave_idx_type& CJ, double* WGHT,const double* WP,
- const octave_idx_type* IWP, double* B, const octave_idx_type& EPLIN,
- octave_idx_type& IER, const double* RPAR, const octave_idx_type* IPAR);
+ const F77_INT& CJ, double* WGHT,const double* WP,
+ const F77_INT* IWP, double* B, const F77_INT& EPLIN,
+ F77_INT& IER, const double* RPAR, const F77_INT* IPAR);
-typedef octave_idx_type (*odepkg_ddaskr_rttype)
- (const octave_idx_type& NEQ, const double& T, const double* Y,
- const double* YP, const octave_idx_type& NRT, const double* RVAL,
- const double* RPAR, const octave_idx_type* IPAR);
+typedef F77_RET_T (*odepkg_ddaskr_rttype)
+ (const F77_INT& NEQ, const double& T, const double* Y,
+ const double* YP, const F77_INT& NRT, const double* RVAL,
+ const double* RPAR, const F77_INT* IPAR);
-// typedef octave_idx_type (*odepkg_ddaskr_krylov_jactype)
-// (odepkg_ddaskr_restype, octave_idx_type& IRES, const octave_idx_type& NEQ,
+// typedef F77_RET_T (*odepkg_ddaskr_krylov_jactype)
+// (odepkg_ddaskr_restype, F77_INT& IRES, const F77_INT& NEQ,
// const double& T, const double* Y, const double* YPRIME,
// double* REWT, const double* SAVR, const double* WK,
// const double& H, const double& CJ, const double* WP,
-// const octave_idx_type* IWP, octave_idx_type& IER,
-// const double* RPAR, const octave_idx_type* IPAR);
+// const F77_INT* IWP, F77_INT& IER,
+// const double* RPAR, const F77_INT* IPAR);
-// typedef octave_idx_type (*odepkg_ddaskr_krylov_psoltype)
-// (const octave_idx_type& NEQ, const double& T, const double* Y,
+// typedef F77_RET_T (*odepkg_ddaskr_krylov_psoltype)
+// (const F77_INT& NEQ, const double& T, const double* Y,
// const double* YPRIME, const double* SAVR, const double* WK,
-// const octave_idx_type& CJ, double* WGHT,const double* WP,
-// const octave_idx_type* IWP, double* B, const octave_idx_type& EPLIN,
-// octave_idx_type& IER, const double* RPAR, const octave_idx_type* IPAR);
+// const F77_INT& CJ, double* WGHT,const double* WP,
+// const F77_INT* IWP, double* B, const F77_INT& EPLIN,
+// F77_INT& IER, const double* RPAR, const F77_INT* IPAR);
extern "C" {
F77_RET_T F77_FUNC (ddaskr, DDASKR)
- (odepkg_ddaskr_restype, const octave_idx_type& NEQ, const double& T,
+ (odepkg_ddaskr_restype, const F77_INT& NEQ, const double& T,
const double* Y, const double* YPRIME, const double& TOUT,
- const octave_idx_type* INFO, const double* RTOL, const double* ATOL,
- const octave_idx_type& IDID, const double* RWORK, const octave_idx_type& LRW,
- const octave_idx_type* IWORK, const octave_idx_type& LIW, const double* RPAR,
- const octave_idx_type* IPAR, odepkg_ddaskr_jactype, odepkg_ddaskr_psoltype,
- odepkg_ddaskr_rttype, const octave_idx_type& NRT, octave_idx_type* JROOT);
+ const F77_INT* INFO, const double* RTOL, const double* ATOL,
+ const F77_INT& IDID, const double* RWORK, const F77_INT& LRW,
+ const F77_INT* IWORK, const F77_INT& LIW, const double* RPAR,
+ const F77_INT* IPAR, odepkg_ddaskr_jactype, odepkg_ddaskr_psoltype,
+ odepkg_ddaskr_rttype, const F77_INT& NRT, F77_INT* JROOT);
}
static octave_value_list vddaskrextarg;
static octave_value vddaskrodefun;
static octave_value vddaskrjacfun;
-static octave_idx_type vddaskrneqn;
+static F77_INT vddaskrneqn;
-octave_idx_type odepkg_ddaskr_resfcn
+F77_RET_T odepkg_ddaskr_resfcn
(const double& T, const double* Y, const double* YPRIME,
- const double& CJ, double* DELTA, octave_idx_type& IRES,
- const double* RPAR, const octave_idx_type* IPAR) {
+ const double& CJ, double* DELTA, F77_INT& IRES,
+ const double* RPAR, const F77_INT* IPAR) {
// Copy the values that come from the Fortran function element wise,
// otherwise Octave will crash if these variables will be freed
ColumnVector A(vddaskrneqn), APRIME(vddaskrneqn);
- for (octave_idx_type vcnt = 0; vcnt < vddaskrneqn; vcnt++) {
+ for (F77_INT vcnt = 0; vcnt < vddaskrneqn; vcnt++) {
A(vcnt) = Y[vcnt]; APRIME(vcnt) = YPRIME[vcnt];
}
@@ -111,28 +111,28 @@ octave_idx_type odepkg_ddaskr_resfcn
// function that keeps the set of implicit differential equations
octave_value_list varin;
varin(0) = T; varin(1) = A; varin(2) = APRIME;
- for (octave_idx_type vcnt = 0; vcnt < vddaskrextarg.length (); vcnt++)
+ for (F77_INT vcnt = 0; vcnt < vddaskrextarg.length (); vcnt++)
varin(vcnt+3) = vddaskrextarg(vcnt);
octave_value_list vout = feval (vddaskrodefun.function_value (), varin, 1);
// Return the results from the function evaluation to the Fortran
// solver, again copy them and don't just create a Fortran vector
ColumnVector vcol = vout(0).column_vector_value ();
- for (octave_idx_type vcnt = 0; vcnt < vddaskrneqn; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < vddaskrneqn; vcnt++)
DELTA[vcnt] = vcol(vcnt);
- return (true);
+ F77_RETURN (true);
}
-octave_idx_type odepkg_ddaskr_jacfcn
+F77_RET_T odepkg_ddaskr_jacfcn
(const double& T, const double* Y, const double* YPRIME,
double* PD, const double& CJ, const double* RPAR,
- const octave_idx_type* IPAR) {
+ const F77_INT* IPAR) {
// Copy the values that come from the Fortran function element-wise,
// otherwise Octave will crash if these variables are freed
ColumnVector A(vddaskrneqn), APRIME(vddaskrneqn);
- for (octave_idx_type vcnt = 0; vcnt < vddaskrneqn; vcnt++) {
+ for (F77_INT vcnt = 0; vcnt < vddaskrneqn; vcnt++) {
A(vcnt) = Y[vcnt]; APRIME(vcnt) = YPRIME[vcnt];
}
@@ -148,37 +148,37 @@ octave_idx_type odepkg_ddaskr_jacfcn
// Fortran solver of the form PD=DG/DY+1/CON*(DG/DY')
octave_value vbdov = vout(0) + CJ * vout(1);
Matrix vbd = vbdov.matrix_value ();
- for (octave_idx_type vrow = 0; vrow < vddaskrneqn; vrow++)
- for (octave_idx_type vcol = 0; vcol < vddaskrneqn; vcol++)
+ for (F77_INT vrow = 0; vrow < vddaskrneqn; vrow++)
+ for (F77_INT vcol = 0; vcol < vddaskrneqn; vcol++)
PD[vrow+vcol*vddaskrneqn] = vbd (vrow, vcol);
// Don't know what my mistake is but the following code line never
// worked for me (ie. the solver crashes Octave if it is used)
// PD = vbd.fortran_vec ();
- return (true);
-}
-
-octave_idx_type odepkg_ddaskr_rtfcn // this is a dummy function
- ( const octave_idx_type& NEQ, const double& T,
- const double* Y, const double* YP,
- const octave_idx_type& NRT, const double* RVAL,
- const double* RPAR, const octave_idx_type* IPAR) {
- return (true);
+ F77_RETURN (true);
}
-octave_idx_type odepkg_ddaskr_psolfcn // this is a dummy function
- ( const octave_idx_type& NEQ, const double& T,
+F77_RET_T odepkg_ddaskr_rtfcn // this is a dummy function
+ ( const F77_INT& NEQ, const double& T,
+ const double* Y, const double* YP,
+ const F77_INT& NRT, const double* RVAL,
+ const double* RPAR, const F77_INT* IPAR) {
+ F77_RETURN (true);
+}
+
+F77_RET_T odepkg_ddaskr_psolfcn // this is a dummy function
+ ( const F77_INT& NEQ, const double& T,
const double* Y, const double* YPRIME,
const double* SAVR, const double* WK,
- const octave_idx_type& CJ, double* WGHT,
- const double* WP, const octave_idx_type* IWP,
- double* B, const octave_idx_type& EPLIN,
- octave_idx_type& IER, const double* RPAR,
- const octave_idx_type* IPAR) {
- return (true);
+ const F77_INT& CJ, double* WGHT,
+ const double* WP, const F77_INT* IWP,
+ double* B, const F77_INT& EPLIN,
+ F77_INT& IER, const double* RPAR,
+ const F77_INT* IPAR) {
+ F77_RETURN (true);
}
-octave_idx_type odepkg_ddaskr_error (octave_idx_type verr) {
+F77_RET_T odepkg_ddaskr_error (F77_INT verr) {
switch (verr)
{
@@ -191,7 +191,7 @@ octave_idx_type odepkg_ddaskr_error (oct
break;
}
- return (true);
+ F77_RETURN (true);
}
// PKG_ADD: autoload ("odekdi", "dldsolver.oct");
@@ -362,7 +362,7 @@ odekdi (@@odepkg_equations_ilorenz, [0,
// Setting the tolerance type that depends on the types (scalar or
// vector) of the options RelTol and AbsTol
- octave_idx_type vitol;
+ F77_INT vitol;
if (vreltol.is_scalar_type () && vabstol.is_scalar_type ()) vitol = 0;
else if (vreltol.length () == vabstol.length ()) vitol = 1;
else {
@@ -523,33 +523,33 @@ odekdi (@@odepkg_equations_ilorenz, [0,
NDArray vRTOL = vreltol.array_value ();
NDArray vATOL = vabstol.array_value ();
- vddaskrneqn = args(2).length ();
+ vddaskrneqn = TO_F77_INT (args(2).length ());
double T = vTIME(0);
double TEND = vTIME(vTIME.numel () - 1);
double *Y = vY.fortran_vec ();
double *YPRIME = vYPRIME.fortran_vec ();
- octave_idx_type IDID = 0;
+ F77_INT IDID = 0;
double *RTOL = vRTOL.fortran_vec ();
double *ATOL = vATOL.fortran_vec ();
double RPAR[1] = {0.0};
- octave_idx_type IPAR[1] = {0};
+ F77_INT IPAR[1] = {0};
- octave_idx_type NRT = 0;
- OCTAVE_LOCAL_BUFFER (octave_idx_type, JROOT, NRT);
- for (octave_idx_type vcnt = 0; vcnt < NRT; vcnt++) JROOT[vcnt] = 0;
+ F77_INT NRT = 0;
+ OCTAVE_LOCAL_BUFFER (F77_INT, JROOT, NRT);
+ for (F77_INT vcnt = 0; vcnt < NRT; vcnt++) JROOT[vcnt] = 0;
- octave_idx_type LRW = 60 + vddaskrneqn * (9 + vddaskrneqn) + 3 * NRT;
+ F77_INT LRW = 60 + vddaskrneqn * (9 + vddaskrneqn) + 3 * NRT;
OCTAVE_LOCAL_BUFFER (double, RWORK, LRW);
- for (octave_idx_type vcnt = 0; vcnt < LRW; vcnt++) RWORK[vcnt] = 0.0;
+ for (F77_INT vcnt = 0; vcnt < LRW; vcnt++) RWORK[vcnt] = 0.0;
- octave_idx_type LIW = 40 + vddaskrneqn;
- OCTAVE_LOCAL_BUFFER (octave_idx_type, IWORK, LIW);
- for (octave_idx_type vcnt = 0; vcnt < LIW; vcnt++) IWORK[vcnt] = 0;
+ F77_INT LIW = 40 + vddaskrneqn;
+ OCTAVE_LOCAL_BUFFER (F77_INT, IWORK, LIW);
+ for (F77_INT vcnt = 0; vcnt < LIW; vcnt++) IWORK[vcnt] = 0;
- octave_idx_type N = 20;
- OCTAVE_LOCAL_BUFFER (octave_idx_type, INFO, N);
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++) INFO[vcnt] = 0;
+ F77_INT N = 20;
+ OCTAVE_LOCAL_BUFFER (F77_INT, INFO, N);
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++) INFO[vcnt] = 0;
INFO[0] = 0; // Define that it is the initial first call
INFO[1] = vitol; // RelTol/AbsTol are scalars or vectors
@@ -618,7 +618,7 @@ odekdi (@@odepkg_equations_ilorenz, [0,
// This call of the Fortran solver has been successful so let us
// plot the output and save the results
- for (octave_idx_type vcnt = 0; vcnt < vddaskrneqn; vcnt++) {
+ for (F77_INT vcnt = 0; vcnt < vddaskrneqn; vcnt++) {
vcres(vcnt) = Y[vcnt]; vydrs(vcnt) = YPRIME[vcnt];
}
vsol = vcres; vyds = vydrs; vtim = T;
@@ -655,7 +655,7 @@ odekdi (@@odepkg_equations_ilorenz, [0,
// Set up values that come from the last Fortran call and that are
// needed to call the OdePkg output function a last time again
- for (octave_idx_type vcnt = 0; vcnt < vddaskrneqn; vcnt++) {
+ for (F77_INT vcnt = 0; vcnt < vddaskrneqn; vcnt++) {
vcres(vcnt) = Y[vcnt]; vydrs(vcnt) = YPRIME[vcnt];
}
vsol = vcres; vyds = vydrs; vtim = T;
diff --git a/src/odepkg_octsolver_mebdfdae.cc b/src/odepkg_octsolver_mebdfdae.cc
--- a/src/odepkg_octsolver_mebdfdae.cc
+++ b/src/odepkg_octsolver_mebdfdae.cc
@@ -46,32 +46,32 @@ similiar to the implementation of this f
#include
#include "odepkg_auxiliary_functions.h"
-typedef octave_idx_type (*odepkg_mebdfdae_usrtype)
- (const octave_idx_type& N, const double& T, const double* Y,
- double* YDOT, const octave_idx_type* IPAR, const double* RPAR,
- const octave_idx_type& IERR);
+typedef F77_RET_T (*odepkg_mebdfdae_usrtype)
+ (const F77_INT& N, const double& T, const double* Y,
+ double* YDOT, const F77_INT* IPAR, const double* RPAR,
+ const F77_INT& IERR);
-typedef octave_idx_type (*odepkg_mebdfdae_jactype)
+typedef F77_RET_T (*odepkg_mebdfdae_jactype)
(const double& T, const double* Y, double* PD,
- const octave_idx_type& N, const octave_idx_type& MEBAND,
- const octave_idx_type* IPAR, const double* RPAR,
- const octave_idx_type& IERR);
+ const F77_INT& N, const F77_INT& MEBAND,
+ const F77_INT* IPAR, const double* RPAR,
+ const F77_INT& IERR);
-typedef octave_idx_type (*odepkg_mebdfdae_masstype)
- (const octave_idx_type& N, double* AM, const octave_idx_type* MASBND,
- const double* RPAR, const octave_idx_type* IPAR, const octave_idx_type& IERR);
+typedef F77_RET_T (*odepkg_mebdfdae_masstype)
+ (const F77_INT& N, double* AM, const F77_INT* MASBND,
+ const double* RPAR, const F77_INT* IPAR, const F77_INT& IERR);
extern "C" {
F77_RET_T F77_FUNC (mebdf, MEBDF) // 3 arguments per line
- (const octave_idx_type& N, const double& T0, const double& HO,
+ (const F77_INT& N, const double& T0, const double& HO,
const double* Y0, const double& TOUT, const double& TEND,
- const octave_idx_type& MF, octave_idx_type& IDID, const octave_idx_type& LOUT,
- const octave_idx_type& LWORK, const double* WORK, const octave_idx_type& LIWORK,
- const octave_idx_type* IWORK, const octave_idx_type* MBND, const octave_idx_type* MASBND,
- const octave_idx_type& MAXDER, const octave_idx_type& ITOL, const double* RTOL,
- const double* ATOL, const double* RPAR, const octave_idx_type* IPAR,
+ const F77_INT& MF, F77_INT& IDID, const F77_INT& LOUT,
+ const F77_INT& LWORK, const double* WORK, const F77_INT& LIWORK,
+ const F77_INT* IWORK, const F77_INT* MBND, const F77_INT* MASBND,
+ const F77_INT& MAXDER, const F77_INT& ITOL, const double* RTOL,
+ const double* ATOL, const double* RPAR, const F77_INT* IPAR,
odepkg_mebdfdae_usrtype, odepkg_mebdfdae_jactype, odepkg_mebdfdae_masstype,
- octave_idx_type& IERR);
+ F77_INT& IERR);
} // extern "C"
static octave_value_list vmebdfdaeextarg;
@@ -80,43 +80,43 @@ static octave_value vmebdfdaejacfun;
static octave_value vmebdfdaemass;
static octave_value vmebdfdaemassstate;
-octave_idx_type odepkg_mebdfdae_usrfcn
- (const octave_idx_type& N, const double& T, const double* Y,
- double* YDOT, const octave_idx_type* IPAR,
- const double* RPAR, const octave_idx_type& IERR) {
+F77_RET_T odepkg_mebdfdae_usrfcn
+ (const F77_INT& N, const double& T, const double* Y,
+ double* YDOT, const F77_INT* IPAR,
+ const double* RPAR, const F77_INT& IERR) {
// Copy the values that come from the Fortran function element wise,
// otherwise Octave will crash if these variables will be freed
ColumnVector A(N);
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++)
A(vcnt) = Y[vcnt];
// Fill the variable for the input arguments before evaluating the
// function that keeps the set of implicit differential equations
octave_value_list varin;
varin(0) = T; varin(1) = A;
- for (octave_idx_type vcnt = 0; vcnt < vmebdfdaeextarg.length (); vcnt++)
+ for (F77_INT vcnt = 0; vcnt < vmebdfdaeextarg.length (); vcnt++)
varin(vcnt+2) = vmebdfdaeextarg(vcnt);
octave_value_list vout = feval (vmebdfdaeodefun.function_value (), varin, 1);
// Return the results from the function evaluation to the Fortran
// solver, again copy them and don't just create a Fortran vector
ColumnVector vcol = vout(0).column_vector_value ();
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++)
YDOT[vcnt] = vcol(vcnt);
- return (true);
+ F77_RETURN (true);
}
-octave_idx_type odepkg_mebdfdae_jacfcn
- (const double& T, const double* Y, double* PD, const octave_idx_type& N,
- const octave_idx_type& MEBAND, const octave_idx_type* IPAR,
- const double* RPAR, const octave_idx_type& IERR) {
+F77_RET_T odepkg_mebdfdae_jacfcn
+ (const double& T, const double* Y, double* PD, const F77_INT& N,
+ const F77_INT& MEBAND, const F77_INT* IPAR,
+ const double* RPAR, const F77_INT& IERR) {
// Copy the values that come from the Fortran function element-wise,
// otherwise Octave will crash if these variables are freed
ColumnVector A(N);
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++)
A(vcnt) = Y[vcnt];
// Set the values that are needed as input arguments before calling
@@ -130,25 +130,25 @@ octave_idx_type odepkg_mebdfdae_jacfcn
// Fortran solver of the form PD=DG/DY+1/CON*(DG/DY')
// octave_value vbdov = vout(0) + 1/CON * vout(1);
Matrix vbd = vout(0).matrix_value ();
- for (octave_idx_type vrow = 0; vrow < N; vrow++)
- for (octave_idx_type vcol = 0; vcol < N; vcol++)
+ for (F77_INT vrow = 0; vrow < N; vrow++)
+ for (F77_INT vcol = 0; vcol < N; vcol++)
PD[vrow+vcol*N] = vbd (vrow, vcol);
// Don't know what my mistake is but the following code line never
// worked for me (ie. the solver crashes Octave if it is used)
// PD = vbd.fortran_vec ();
- return (true);
+ F77_RETURN (true);
}
-octave_idx_type odepkg_mebdfdae_massfcn
- (const octave_idx_type& N, double* AM,
- const octave_idx_type* MASBND, const double* RPAR,
- const octave_idx_type* IPAR, const octave_idx_type& IERR) {
+F77_RET_T odepkg_mebdfdae_massfcn
+ (const F77_INT& N, double* AM,
+ const F77_INT* MASBND, const double* RPAR,
+ const F77_INT* IPAR, const F77_INT& IERR) {
// Copy the values that come from the Fortran function element-wise,
// otherwise Octave will crash if these variables are freed
ColumnVector A(N);
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++)
A(vcnt) = 0.0;
// Set the values that are needed as input arguments before calling
@@ -159,14 +159,14 @@ octave_idx_type odepkg_mebdfdae_massfcn
(vmebdfdaemass, vmebdfdaemassstate, vt, vy, vmebdfdaeextarg);
Matrix vam = vout.matrix_value ();
- for (octave_idx_type vrow = 0; vrow < N; vrow++)
- for (octave_idx_type vcol = 0; vcol < N; vcol++) {
+ for (F77_INT vrow = 0; vrow < N; vrow++)
+ for (F77_INT vcol = 0; vcol < N; vcol++) {
AM[vrow+vcol*N] = vam (vrow, vcol);
}
- return (true);
+ F77_RETURN (true);
}
-octave_idx_type odepkg_mebdfdae_error (octave_idx_type verr) {
+F77_RET_T odepkg_mebdfdae_error (F77_INT verr) {
switch (verr)
{
@@ -238,7 +238,7 @@ error number \"%d\")", verr);
break;
}
- return (true);
+ F77_RETURN (true);
}
// PKG_ADD: autoload ("odebda", "dldsolver.oct");
@@ -394,7 +394,7 @@ odebda (@@odepkg_equations_lorenz, [0, 2
// Setting the tolerance type that depends on the types (scalar or
// vector) of the options RelTol and AbsTol
- octave_idx_type vitol = 0;
+ F77_INT vitol = 0;
if (vreltol.is_scalar_type () && vabstol.is_scalar_type ()) vitol = 2;
else if (vreltol.is_scalar_type () && !vabstol.is_scalar_type ()) vitol = 3;
else if (!vreltol.is_scalar_type () && vabstol.is_scalar_type ()) vitol = 4;
@@ -469,7 +469,7 @@ odebda (@@odepkg_equations_lorenz, [0, 2
// The options 'Jacobian', 'JPattern' and 'Vectorized'
octave_value vjac = vodeopt.contents ("Jacobian");
- octave_idx_type vmebdfdaejac = 22; // We need to set this if no Jac available
+ F77_INT vmebdfdaejac = 22; // We need to set this if no Jac available
if (!vjac.is_empty ()) {
vmebdfdaejacfun = vjac; vmebdfdaejac = 21;
}
@@ -477,7 +477,7 @@ odebda (@@odepkg_equations_lorenz, [0, 2
// Implementation of the option 'Mass' has been finished, these
// options can be set by the user to another value than default
vmebdfdaemass = vodeopt.contents ("Mass");
- octave_idx_type vmebdfdaemas = 0;
+ F77_INT vmebdfdaemas = 0;
if (!vmebdfdaemass.is_empty ()) {
vmebdfdaemas = 1;
if (vmebdfdaemass.is_function_handle () || vmebdfdaemass.is_inline_function ())
@@ -560,34 +560,34 @@ odebda (@@odepkg_equations_lorenz, [0, 2
NDArray vATOL = vabstol.array_value ();
NDArray vY0 = args(2).array_value ();
- octave_idx_type N = args(2).length ();
+ F77_INT N = TO_F77_INT (args(2).length ());
double T0 = vTIME(0);
double HO = vinitstep.double_value ();
double *Y0 = vY0.fortran_vec ();
double TOUT = T0 + vinitstep.double_value ();
double TEND = vTIME(vTIME.numel ()-1);
- octave_idx_type MF = vmebdfdaejac;
- octave_idx_type IDID = 1;
- octave_idx_type LOUT = 6; // Logical output channel "not opened"
- octave_idx_type LWORK = 42*N+3*N*N+4;
+ F77_INT MF = vmebdfdaejac;
+ F77_INT IDID = 1;
+ F77_INT LOUT = 6; // Logical output channel "not opened"
+ F77_INT LWORK = 42*N+3*N*N+4;
OCTAVE_LOCAL_BUFFER (double, WORK, LWORK);
- for (octave_idx_type vcnt = 0; vcnt < LWORK; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < LWORK; vcnt++)
WORK[vcnt] = 0.0;
- octave_idx_type LIWORK = N+14;
- OCTAVE_LOCAL_BUFFER (octave_idx_type, IWORK, LIWORK);
- for (octave_idx_type vcnt = 0; vcnt < LIWORK; vcnt++)
+ F77_INT LIWORK = N+14;
+ OCTAVE_LOCAL_BUFFER (F77_INT, IWORK, LIWORK);
+ for (F77_INT vcnt = 0; vcnt < LIWORK; vcnt++)
IWORK[vcnt] = 0;
- octave_idx_type MBND[4] = {N, N, N, N};
- octave_idx_type MASBND[4] = {0, N, 0, N};
+ F77_INT MBND[4] = {N, N, N, N};
+ F77_INT MASBND[4] = {0, N, 0, N};
if (!vmebdfdaemass.is_empty ()) MASBND[0] = 1;
- octave_idx_type MAXDER = vmaxder.int_value ();
- octave_idx_type ITOL = vitol;
+ F77_INT MAXDER = vmaxder.int_value ();
+ F77_INT ITOL = vitol;
double *RTOL = vRTOL.fortran_vec ();
double *ATOL = vATOL.fortran_vec ();
double RPAR[1] = {0.0};
- octave_idx_type IPAR[1] = {0};
- octave_idx_type IERR = 0;
+ F77_INT IPAR[1] = {0};
+ F77_INT IERR = 0;
IWORK[0] = N; // Number of variables of index 1
IWORK[1] = 0; // Number of variables of index 2
@@ -648,7 +648,7 @@ odebda (@@odepkg_equations_lorenz, [0, 2
// This call of the Fortran solver has been successful so let us
// plot the output and save the results
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++)
vcres(vcnt) = Y0[vcnt];
vsol = vcres; vtim = TOUT;
if (!vevents.is_empty ()) {
@@ -672,8 +672,8 @@ odebda (@@odepkg_equations_lorenz, [0, 2
}
}
else { // if (vTIME.numel () > 2) we have all the time values needed
- volatile octave_idx_type vtimecnt = 1;
- octave_idx_type vtimelen = vTIME.numel ();
+ volatile F77_INT vtimecnt = 1;
+ F77_INT vtimelen = TO_F77_INT (vTIME.numel ());
while (vtimecnt < vtimelen) {
vtimecnt++; TOUT = vTIME(vtimecnt-1);
@@ -695,7 +695,7 @@ odebda (@@odepkg_equations_lorenz, [0, 2
// The last call of the Fortran solver has been successful so
// let us plot the output and save the results
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++)
vcres(vcnt) = Y0[vcnt];
vsol = vcres; vtim = TOUT;
if (!vevents.is_empty ()) {
@@ -725,7 +725,7 @@ odebda (@@odepkg_equations_lorenz, [0, 2
// Set up values that come from the last Fortran call and that are
// needed to call the OdePkg output function one last time again
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++) vcres(vcnt) = Y0[vcnt];
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++) vcres(vcnt) = Y0[vcnt];
vsol = vcres; vtim = TOUT;
if (!vevents.is_empty ())
odepkg_auxiliary_evaleventfun (vevents, vtim, vsol, vmebdfdaeextarg, 2);
diff --git a/src/odepkg_octsolver_mebdfi.cc b/src/odepkg_octsolver_mebdfi.cc
--- a/src/odepkg_octsolver_mebdfi.cc
+++ b/src/odepkg_octsolver_mebdfi.cc
@@ -45,59 +45,59 @@ will be broken, eg.
/* -*- texinfo -*-
* @subsection Source File @file{odepkg_octsolver_mebdfi.cc}
*
- * @deftp {Typedef} {octave_idx_type (*odepkg_mebdfi_usrtype)}
+ * @deftp {Typedef} {F77_INT (*odepkg_mebdfi_usrtype)}
* This @code{typedef} is used to define the input and output arguments of the user function for the IDE problem that is further needed by the Fortran core solver @code{mebdfi}. The implementation of this @code{typedef} is
*
* @example
- * typedef octave_idx_type (*odepkg_mebdfi_usrtype)
- * (const octave_idx_type& N, const double& T, const double* Y,
- * double* DELTA, const double* YPRIME, const octave_idx_type* IPAR,
- * const double* RPAR, const octave_idx_type& IERR);
+ * typedef F77_RET_T (*odepkg_mebdfi_usrtype)
+ * (const F77_INT& N, const double& T, const double* Y,
+ * double* DELTA, const double* YPRIME, const F77_INT* IPAR,
+ * const double* RPAR, const F77_INT& IERR);
* @end example
* @end deftp
*/
-typedef octave_idx_type (*odepkg_mebdfi_usrtype)
- (const octave_idx_type& N, const double& T, const double* Y,
- double* DELTA, const double* YPRIME, const octave_idx_type* IPAR,
- const double* RPAR, const octave_idx_type& IERR);
+typedef F77_RET_T (*odepkg_mebdfi_usrtype)
+ (const F77_INT& N, const double& T, const double* Y,
+ double* DELTA, const double* YPRIME, const F77_INT* IPAR,
+ const double* RPAR, const F77_INT& IERR);
/* -*- texinfo -*-
- * @deftp {Typedef} {octave_idx_type (*odepkg_mebdfi_jactype)}
+ * @deftp {Typedef} {F77_INT (*odepkg_mebdfi_jactype)}
*
* This @code{typedef} is used to define the input and output arguments of the @code{Jacobian} function for the IDE problem that is further needed by the Fortran core solver @code{mebdfi}. The implementation of this @code{typedef} is
*
* @example
- * typedef octave_idx_type (*odepkg_mebdfi_jactype)
+ * typedef F77_RET_T (*odepkg_mebdfi_jactype)
* (const double& T, const double* Y, double* PD,
- * const octave_idx_type& N, const double* YPRIME,
- * const octave_idx_type* MBND, const double& CON,
- * const octave_idx_type* IPAR, const double* RPAR,
- * const octave_idx_type& IERR);
+ * const F77_INT& N, const double* YPRIME,
+ * const F77_INT* MBND, const double& CON,
+ * const F77_INT* IPAR, const double* RPAR,
+ * const F77_INT& IERR);
* @end example
* @end deftp
*/
-typedef octave_idx_type (*odepkg_mebdfi_jactype) // 3 arguments per line
+typedef F77_RET_T (*odepkg_mebdfi_jactype) // 3 arguments per line
(const double& T, const double* Y, double* PD,
- const octave_idx_type& N, const double* YPRIME, const octave_idx_type* MBND,
- const double& CON, const octave_idx_type* IPAR, const double* RPAR,
- const octave_idx_type& IERR);
+ const F77_INT& N, const double* YPRIME, const F77_INT* MBND,
+ const double& CON, const F77_INT* IPAR, const double* RPAR,
+ const F77_INT& IERR);
extern "C" {
/* -*- texinfo -*-
- * @deftp {Prototype} {F77_RET_T F77_FUNC (mebdfi, MEBDFI)} (const octave_idx_type& N, const double& T0, const double& HO, const double* Y0, const double* YPRIME, const double& TOUT, const double& TEND, const octave_idx_type& MF, octave_idx_type& IDID, const octave_idx_type& LOUT, const octave_idx_type& LWORK, const double* WORK, const octave_idx_type& LIWORK, const octave_idx_type* IWORK, const octave_idx_type* MBND, const octave_idx_type& MAXDER, const octave_idx_type& ITOL, const double* RTOL, const double* ATOL, const double* RPAR, const octave_idx_type* IPAR, odepkg_mebdfi_jactype, odepkg_mebdfi_usrtype, octave_idx_type& IERR);
+ * @deftp {Prototype} {F77_RET_T F77_FUNC (mebdfi, MEBDFI)} (const F77_INT& N, const double& T0, const double& HO, const double* Y0, const double* YPRIME, const double& TOUT, const double& TEND, const F77_INT& MF, F77_INT& IDID, const F77_INT& LOUT, const F77_INT& LWORK, const double* WORK, const F77_INT& LIWORK, const F77_INT* IWORK, const F77_INT* MBND, const F77_INT& MAXDER, const F77_INT& ITOL, const double* RTOL, const double* ATOL, const double* RPAR, const F77_INT* IPAR, odepkg_mebdfi_jactype, odepkg_mebdfi_usrtype, F77_INT& IERR);
*
* The prototype @code{F77_FUNC (mebdfi, MEBDFI)} is used to represent the information about the Fortran core solver @code{mebdfi} that is defined in the Fortran source file @file{mebdfi.f} (cf. the Fortran source file @file{mebdfi.f} for further details).
* @end deftp
*/
F77_RET_T F77_FUNC (mebdfi, MEBDFI) // 3 arguments per line
- (const octave_idx_type& N, const double& T0, const double& HO,
+ (const F77_INT& N, const double& T0, const double& HO,
const double* Y0, const double* YPRIME, const double& TOUT,
- const double& TEND, const octave_idx_type& MF, octave_idx_type& IDID,
- const octave_idx_type& LOUT, const octave_idx_type& LWORK, const double* WORK,
- const octave_idx_type& LIWORK, const octave_idx_type* IWORK, const octave_idx_type* MBND,
- const octave_idx_type& MAXDER, const octave_idx_type& ITOL, const double* RTOL,
- const double* ATOL, const double* RPAR, const octave_idx_type* IPAR,
- odepkg_mebdfi_jactype, odepkg_mebdfi_usrtype, octave_idx_type& IERR);
+ const double& TEND, const F77_INT& MF, F77_INT& IDID,
+ const F77_INT& LOUT, const F77_INT& LWORK, const double* WORK,
+ const F77_INT& LIWORK, const F77_INT* IWORK, const F77_INT* MBND,
+ const F77_INT& MAXDER, const F77_INT& ITOL, const double* RTOL,
+ const double* ATOL, const double* RPAR, const F77_INT* IPAR,
+ odepkg_mebdfi_jactype, odepkg_mebdfi_usrtype, F77_INT& IERR);
}
/* -*- texinfo -*-
@@ -125,7 +125,7 @@ static octave_value vmebdfiodefun;
static octave_value vmebdfijacfun;
/* -*- texinfo -*-
- * @deftypefn {Function} {octave_idx_type} {odepkg_mebdfi_usrfcn} (const octave_idx_type& N, const double& T, const double* Y, double* DELTA, const double* YPRIME, const octave_idx_type* IPAR, const double* RPAR, const octave_idx_type& IERR)
+ * @deftypefn {Function} {F77_INT} {odepkg_mebdfi_usrfcn} (const F77_INT& N, const double& T, const double* Y, double* DELTA, const double* YPRIME, const F77_INT* IPAR, const double* RPAR, const F77_INT& IERR)
*
* Return @code{true} if the evaluation of the user function was successful, return @code{false} otherwise. This function is directly called from the Fortran core solver @code{mebdfi}. The input arguments of this function are
*
@@ -141,15 +141,15 @@ static octave_value vmebdfijacfun;
* @end itemize
* @end deftypefn
*/
-octave_idx_type odepkg_mebdfi_usrfcn
- (const octave_idx_type& N, const double& T, const double* Y,
- double* DELTA, const double* YPRIME, const octave_idx_type* IPAR,
- const double* RPAR, const octave_idx_type& IERR) {
+F77_RET_T odepkg_mebdfi_usrfcn
+ (const F77_INT& N, const double& T, const double* Y,
+ double* DELTA, const double* YPRIME, const F77_INT* IPAR,
+ const double* RPAR, const F77_INT& IERR) {
// Copy the values that come from the Fortran function element wise,
// otherwise Octave will crash if these variables will be freed
ColumnVector A(N), APRIME(N);
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++) {
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++) {
A(vcnt) = Y[vcnt]; APRIME(vcnt) = YPRIME[vcnt];
}
@@ -157,21 +157,21 @@ octave_idx_type odepkg_mebdfi_usrfcn
// function that keeps the set of implicit differential equations
octave_value_list varin;
varin(0) = T; varin(1) = A; varin(2) = APRIME;
- for (octave_idx_type vcnt = 0; vcnt < vmebdfiextarg.length (); vcnt++)
+ for (F77_INT vcnt = 0; vcnt < vmebdfiextarg.length (); vcnt++)
varin(vcnt+3) = vmebdfiextarg(vcnt);
octave_value_list vout = feval (vmebdfiodefun.function_value (), varin, 1);
// Return the results from the function evaluation to the Fortran
// solver, again copy them and don't just create a Fortran vector
ColumnVector vcol = vout(0).column_vector_value ();
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++)
DELTA[vcnt] = vcol(vcnt);
- return (true);
+ F77_RETURN (true);
}
/* -*- texinfo -*-
- * @deftypefn {Function} {octave_idx_type} {odepkg_mebdfi_jacfcn} (const double& T, const double* Y, double* PD, const octave_idx_type& N, const double* YPRIME, const octave_idx_type* MBND, const double& CON, const octave_idx_type* IPAR, const double* RPAR, const octave_idx_type& IERR)
+ * @deftypefn {Function} {F77_INT} {odepkg_mebdfi_jacfcn} (const double& T, const double* Y, double* PD, const F77_INT& N, const double* YPRIME, const F77_INT* MBND, const double& CON, const F77_INT* IPAR, const double* RPAR, const F77_INT& IERR)
*
* Return @code{true} if the evaluation of the Jacobian function (that is defined for a special IDE problem in Octave) was successful, otherwise return @code{false}. This function is directly called from the Fortran core solver @code{mebdfi}. The input arguments of this function are
*
@@ -189,16 +189,16 @@ octave_idx_type odepkg_mebdfi_usrfcn
* @end itemize
* @end deftypefn
*/
-octave_idx_type odepkg_mebdfi_jacfcn
- (const double& T, const double* Y, double* PD, const octave_idx_type& N,
- const double* YPRIME, const octave_idx_type* MBND,
- const double& CON, const octave_idx_type* IPAR,
- const double* RPAR, const octave_idx_type& IERR) {
+F77_RET_T odepkg_mebdfi_jacfcn
+ (const double& T, const double* Y, double* PD, const F77_INT& N,
+ const double* YPRIME, const F77_INT* MBND,
+ const double& CON, const F77_INT* IPAR,
+ const double* RPAR, const F77_INT& IERR) {
// Copy the values that come from the Fortran function element-wise,
// otherwise Octave will crash if these variables are freed
ColumnVector A(N), APRIME(N);
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++) {
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++) {
A(vcnt) = Y[vcnt]; APRIME(vcnt) = YPRIME[vcnt];
}
@@ -214,22 +214,22 @@ octave_idx_type odepkg_mebdfi_jacfcn
// Fortran solver of the form PD=DG/DY+1/CON*(DG/DY')
octave_value vbdov = vout(0) + 1/CON * vout(1);
Matrix vbd = vbdov.matrix_value ();
- for (octave_idx_type vrow = 0; vrow < N; vrow++)
- for (octave_idx_type vcol = 0; vcol < N; vcol++)
+ for (F77_INT vrow = 0; vrow < N; vrow++)
+ for (F77_INT vcol = 0; vcol < N; vcol++)
PD[vrow+vcol*N] = vbd (vrow, vcol);
// Don't know what my mistake is but the following code line never
// worked for me (ie. the solver crashes Octave if it is used)
// PD = vbd.fortran_vec ();
- return (true);
+ F77_RETURN (true);
}
/* -*- texinfo -*-
- * @deftypefn {Function} octave_idx_type odepkg_mebdfi_error (octave_idx_type verr)
+ * @deftypefn {Function} F77_INT odepkg_mebdfi_error (F77_INT verr)
* TODO
* @end deftypefn
*/
-octave_idx_type odepkg_mebdfi_error (octave_idx_type verr) {
+F77_RET_T odepkg_mebdfi_error (F77_INT verr) {
switch (verr)
{
@@ -301,7 +301,7 @@ error number \"%d\")", verr);
break;
}
- return (true);
+ F77_RETURN (true);
}
/* -*- texinfo -*-
@@ -486,7 +486,7 @@ odebdi (@@odepkg_equations_ilorenz, [0,
// Setting the tolerance type that depends on the types (scalar or
// vector) of the options RelTol and AbsTol
- octave_idx_type vitol = 0;
+ F77_INT vitol = 0;
if (vreltol.is_scalar_type () && vabstol.is_scalar_type ()) vitol = 2;
else if (vreltol.is_scalar_type () && !vabstol.is_scalar_type ()) vitol = 3;
else if (!vreltol.is_scalar_type () && vabstol.is_scalar_type ()) vitol = 4;
@@ -561,7 +561,7 @@ odebdi (@@odepkg_equations_ilorenz, [0,
// The options 'Jacobian', 'JPattern' and 'Vectorized'
octave_value vjac = vodeopt.contents ("Jacobian");
- octave_idx_type vmebdfijac = 22; // We need to set this if no Jac available
+ F77_INT vmebdfijac = 22; // We need to set this if no Jac available
if (!vjac.is_empty ()) {
vmebdfijacfun = vjac; vmebdfijac = 21;
}
@@ -649,7 +649,7 @@ odebdi (@@odepkg_equations_ilorenz, [0,
NDArray vY0 = args(2).array_value ();
NDArray vYPRIME = args(3).array_value ();
- octave_idx_type N = args(2).length ();
+ F77_INT N = TO_F77_INT (args(2).length ());
double T0 = vTIME(0);
double HO = vinitstep.double_value ();
double *Y0 = vY0.fortran_vec ();
@@ -657,26 +657,26 @@ odebdi (@@odepkg_equations_ilorenz, [0,
double TOUT = T0 + vinitstep.double_value ();
double TEND = vTIME(vTIME.numel ()-1);
- octave_idx_type MF = vmebdfijac;
- octave_idx_type IDID = 1;
- octave_idx_type LOUT = 6; // Logical output channel "not opened"
- octave_idx_type LWORK = 32*N+2*N*N+3;
+ F77_INT MF = vmebdfijac;
+ F77_INT IDID = 1;
+ F77_INT LOUT = 6; // Logical output channel "not opened"
+ F77_INT LWORK = 32*N+2*N*N+3;
OCTAVE_LOCAL_BUFFER (double, WORK, LWORK);
- for (octave_idx_type vcnt = 0; vcnt < LWORK; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < LWORK; vcnt++)
WORK[vcnt] = 0.0;
- octave_idx_type LIWORK = N+14;
- OCTAVE_LOCAL_BUFFER (octave_idx_type, IWORK, LIWORK);
- for (octave_idx_type vcnt = 0; vcnt < LIWORK; vcnt++)
+ F77_INT LIWORK = N+14;
+ OCTAVE_LOCAL_BUFFER (F77_INT, IWORK, LIWORK);
+ for (F77_INT vcnt = 0; vcnt < LIWORK; vcnt++)
IWORK[vcnt] = 0;
- octave_idx_type MBND[4] = {N, N, N, N};
- octave_idx_type MAXDER = vmaxder.int_value ();
+ F77_INT MBND[4] = {N, N, N, N};
+ F77_INT MAXDER = vmaxder.int_value ();
- octave_idx_type ITOL = vitol;
+ F77_INT ITOL = vitol;
double *RTOL = vRTOL.fortran_vec ();
double *ATOL = vATOL.fortran_vec ();
double RPAR[1] = {0.0};
- octave_idx_type IPAR[1] = {0};
- octave_idx_type IERR = 0;
+ F77_INT IPAR[1] = {0};
+ F77_INT IERR = 0;
IWORK[13] = 1000000; // The maximum number of steps allowed
@@ -740,7 +740,7 @@ odebdi (@@odepkg_equations_ilorenz, [0,
// This call of the Fortran solver has been successful so let us
// plot the output and save the results
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++) {
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++) {
vcres(vcnt) = Y0[vcnt]; vydrs(vcnt) = YPRIME[vcnt];
}
vsol = vcres; vyds = vydrs; vtim = TOUT;
@@ -768,8 +768,8 @@ odebdi (@@odepkg_equations_ilorenz, [0,
}
}
else { // if (vTIME.numel () > 2) we have all the time values needed
- volatile octave_idx_type vtimecnt = 1;
- octave_idx_type vtimelen = vTIME.numel () - 1;
+ volatile F77_INT vtimecnt = 1;
+ F77_INT vtimelen = TO_F77_INT (vTIME.numel () - 1);
while (vtimecnt < vtimelen) {
vtimecnt++; TOUT = vTIME(vtimecnt);
@@ -791,7 +791,7 @@ odebdi (@@odepkg_equations_ilorenz, [0,
// The last call of the Fortran solver has been successful so
// let us plot the output and save the results
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++) {
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++) {
vcres(vcnt) = Y0[vcnt]; vydrs(vcnt) = YPRIME[vcnt];
}
vsol = vcres; vyds = vydrs; vtim = TOUT;
@@ -825,7 +825,7 @@ odebdi (@@odepkg_equations_ilorenz, [0,
// Set up values that come from the last Fortran call and that are
// needed to call the OdePkg output function one last time again
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++) {
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++) {
vcres(vcnt) = Y0[vcnt]; vydrs(vcnt) = YPRIME[vcnt];
}
vsol = vcres; vyds = vydrs; vtim = TOUT;
diff --git a/src/odepkg_octsolver_radau.cc b/src/odepkg_octsolver_radau.cc
--- a/src/odepkg_octsolver_radau.cc
+++ b/src/odepkg_octsolver_radau.cc
@@ -37,46 +37,46 @@ similiar to the implementation of this f
#include
#include "odepkg_auxiliary_functions.h"
-typedef octave_idx_type (*odepkg_radau_usrtype)
- (const octave_idx_type& N, const double& X, const double* Y, double* F,
+typedef F77_RET_T (*odepkg_radau_usrtype)
+ (const F77_INT& N, const double& X, const double* Y, double* F,
const double* RPAR,
- const octave_idx_type* IPAR);
+ const F77_INT* IPAR);
-typedef octave_idx_type (*odepkg_radau_jactype)
- (const octave_idx_type& N, const double& X, const double* Y, double* DFY,
- const octave_idx_type& LDFY,
+typedef F77_RET_T (*odepkg_radau_jactype)
+ (const F77_INT& N, const double& X, const double* Y, double* DFY,
+ const F77_INT& LDFY,
const double* RPAR,
- const octave_idx_type* IPAR);
+ const F77_INT* IPAR);
typedef F77_RET_T (*odepkg_radau_masstype)
- (const octave_idx_type& N, double* AM,
- const octave_idx_type* LMAS,
+ (const F77_INT& N, double* AM,
+ const F77_INT* LMAS,
const double* RPAR,
- const octave_idx_type* IPAR);
+ const F77_INT* IPAR);
-typedef octave_idx_type (*odepkg_radau_soltype)
- (const octave_idx_type& NR, const double& XOLD, const double& X,
- const double* Y, const double* CONT, const octave_idx_type* LRC,
- const octave_idx_type& N, const double* RPAR,
- const octave_idx_type* IPAR, octave_idx_type& IRTRN);
+typedef F77_RET_T (*odepkg_radau_soltype)
+ (const F77_INT& NR, const double& XOLD, const double& X,
+ const double* Y, const double* CONT, const F77_INT* LRC,
+ const F77_INT& N, const double* RPAR,
+ const F77_INT* IPAR, F77_INT& IRTRN);
extern "C" {
F77_RET_T F77_FUNC (radau, RADAU)
- (const octave_idx_type& N, odepkg_radau_usrtype,
+ (const F77_INT& N, odepkg_radau_usrtype,
const double& X, const double* Y, const double& XEND,
const double& H, const double* RTOL, const double* ATOL,
- const octave_idx_type& ITOL, odepkg_radau_jactype, const octave_idx_type& IJAC,
- const octave_idx_type& MLJAC, const octave_idx_type& MUJAC, odepkg_radau_masstype,
- const octave_idx_type& IMAS, const octave_idx_type& MLMAS,
- const octave_idx_type& MUMAS, odepkg_radau_soltype, const octave_idx_type& IOUT,
- const double* WORK, const octave_idx_type& LWORK, const octave_idx_type* IWORK,
- const octave_idx_type& LIWORK, const double* RPAR,
- const octave_idx_type* IPAR, const octave_idx_type& IDID);
+ const F77_INT& ITOL, odepkg_radau_jactype, const F77_INT& IJAC,
+ const F77_INT& MLJAC, const F77_INT& MUJAC, odepkg_radau_masstype,
+ const F77_INT& IMAS, const F77_INT& MLMAS,
+ const F77_INT& MUMAS, odepkg_radau_soltype, const F77_INT& IOUT,
+ const double* WORK, const F77_INT& LWORK, const F77_INT* IWORK,
+ const F77_INT& LIWORK, const double* RPAR,
+ const F77_INT* IPAR, const F77_INT& IDID);
double F77_FUNC (contra, CONTRA)
- (const octave_idx_type& I, const double& S,
- const double* CONT, const octave_idx_type* LRC);
+ (const F77_INT& I, const double& S,
+ const double* CONT, const F77_INT* LRC);
} // extern "C"
@@ -91,15 +91,15 @@ static octave_value vradaurefine;
static octave_value vradaumass;
static octave_value vradaumassstate;
-octave_idx_type odepkg_radau_usrfcn
- (const octave_idx_type& N, const double& X, const double* Y,
+F77_RET_T odepkg_radau_usrfcn
+ (const F77_INT& N, const double& X, const double* Y,
double* F, const double* RPAR,
- const octave_idx_type* IPAR) {
+ const F77_INT* IPAR) {
// Copy the values that come from the Fortran function element wise,
// otherwise Octave will crash if these variables will be freed
ColumnVector A(N);
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++) {
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++) {
A(vcnt) = Y[vcnt];
// octave_stdout << "I am here Y[" << vcnt << "] " << Y[vcnt] << std::endl;
// octave_stdout << "I am here T " << X << std::endl;
@@ -109,29 +109,29 @@ octave_idx_type odepkg_radau_usrfcn
// function that keeps the set of differential algebraic equations
octave_value_list varin;
varin(0) = X; varin(1) = A;
- for (octave_idx_type vcnt = 0; vcnt < vradauextarg.length (); vcnt++)
+ for (F77_INT vcnt = 0; vcnt < vradauextarg.length (); vcnt++)
varin(vcnt+2) = vradauextarg(vcnt);
octave_value_list vout = feval (vradauodefun.function_value (), varin, 1);
// Return the results from the function evaluation to the Fortran
// solver, again copy them and don't just create a Fortran vector
ColumnVector vcol = vout(0).column_vector_value ();
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++)
F[vcnt] = vcol(vcnt);
- return (true);
+ F77_RETURN (true);
}
-octave_idx_type odepkg_radau_jacfcn
- (const octave_idx_type& N, const double& X, const double* Y,
- double* DFY, const octave_idx_type& LDFY,
+F77_RET_T odepkg_radau_jacfcn
+ (const F77_INT& N, const double& X, const double* Y,
+ double* DFY, const F77_INT& LDFY,
const double* RPAR,
- const octave_idx_type* IPAR) {
+ const F77_INT* IPAR) {
// Copy the values that come from the Fortran function element-wise,
// otherwise Octave will crash if these variables are freed
ColumnVector A(N);
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++)
A(vcnt) = Y[vcnt];
// Set the values that are needed as input arguments before calling
@@ -142,23 +142,23 @@ octave_idx_type odepkg_radau_jacfcn
(vradaujacfun, vt, vy, vradauextarg);
Matrix vdfy = vout.matrix_value ();
- for (octave_idx_type vcol = 0; vcol < N; vcol++)
- for (octave_idx_type vrow = 0; vrow < N; vrow++)
+ for (F77_INT vcol = 0; vcol < N; vcol++)
+ for (F77_INT vrow = 0; vrow < N; vrow++)
DFY[vrow+vcol*N] = vdfy (vrow, vcol);
- return (true);
+ F77_RETURN (true);
}
F77_RET_T odepkg_radau_massfcn
- (const octave_idx_type& N, double* AM,
- const octave_idx_type* LMAS,
+ (const F77_INT& N, double* AM,
+ const F77_INT* LMAS,
const double* RPAR,
- const octave_idx_type* IPAR) {
+ const F77_INT* IPAR) {
// Copy the values that come from the Fortran function element-wise,
// otherwise Octave will crash if these variables are freed
ColumnVector A(N);
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++)
A(vcnt) = 0.0;
// warning_with_id ("OdePkg:InvalidFunctionCall",
@@ -172,23 +172,23 @@ F77_RET_T odepkg_radau_massfcn
(vradaumass, vradaumassstate, vt, vy, vradauextarg);
Matrix vam = vout.matrix_value ();
- for (octave_idx_type vrow = 0; vrow < N; vrow++)
- for (octave_idx_type vcol = 0; vcol < N; vcol++)
+ for (F77_INT vrow = 0; vrow < N; vrow++)
+ for (F77_INT vcol = 0; vcol < N; vcol++)
AM[vrow+vcol*N] = vam (vrow, vcol);
return ((F77_RET_T) true);
}
-octave_idx_type odepkg_radau_solfcn
- (const octave_idx_type& NR, const double& XOLD, const double& X,
- const double* Y, const double* CONT, const octave_idx_type* LRC,
- const octave_idx_type& N, const double* RPAR,
- const octave_idx_type* IPAR, octave_idx_type& IRTRN) {
+F77_RET_T odepkg_radau_solfcn
+ (const F77_INT& NR, const double& XOLD, const double& X,
+ const double* Y, const double* CONT, const F77_INT* LRC,
+ const F77_INT& N, const double* RPAR,
+ const F77_INT* IPAR, F77_INT& IRTRN) {
// Copy the values that come from the Fortran function element-wise,
// otherwise Octave will crash if these variables are freed
ColumnVector A(N);
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++)
A(vcnt) = Y[vcnt];
// Set the values that are needed as input arguments before calling
@@ -219,12 +219,12 @@ octave_idx_type odepkg_radau_solfcn
if (!vradaupltfun.is_empty ()) {
if (vradaurefine.int_value () > 0) {
ColumnVector B(N); double vtb = 0.0;
- for (octave_idx_type vcnt = 1; vcnt < vradaurefine.int_value (); vcnt++) {
+ for (F77_INT vcnt = 1; vcnt < vradaurefine.int_value (); vcnt++) {
// Calculate time stamps between XOLD and X and get the
// results at these time stamps
vtb = (X - XOLD) * vcnt / vradaurefine.int_value () + XOLD;
- for (octave_idx_type vcou = 0; vcou < N; vcou++)
+ for (F77_INT vcou = 0; vcou < N; vcou++)
B(vcou) = F77_FUNC (contra, CONTRA) (vcou+1, vtb, CONT, LRC);
// Evaluate the 'OutputFcn' with the approximated values from
@@ -237,11 +237,11 @@ octave_idx_type odepkg_radau_solfcn
}
// Evaluate the 'OutputFcn' with the results from the solver, if
// the OutputFcn returns true then set a negative value in IRTRN
- IRTRN = - odepkg_auxiliary_evalplotfun
- (vradaupltfun, vradauoutsel, vt, vy, vradauextarg, 1);
+ IRTRN = TO_F77_INT (- odepkg_auxiliary_evalplotfun
+ (vradaupltfun, vradauoutsel, vt, vy, vradauextarg, 1));
}
- return (true);
+ F77_RETURN (true);
}
// PKG_ADD: autoload ("ode2r", "dldsolver.oct");
@@ -397,7 +397,7 @@ ode2r (@@odepkg_equations_lorenz, [0, 25
// Setting the tolerance type that depends on the types (scalar or
// vector) of the options RelTol and AbsTol
- octave_idx_type vitol = 0;
+ F77_INT vitol = 0;
if (vreltol.is_scalar_type () && (vreltol.length () == vabstol.length ()))
vitol = 0;
else if (!vreltol.is_scalar_type () && (vreltol.length () == vabstol.length ()))
@@ -472,7 +472,7 @@ ode2r (@@odepkg_equations_lorenz, [0, 25
// Implementation of the option 'Jacobian' has been finished, these
// options can be set by the user to another value than default
vradaujacfun = vodeopt.contents ("Jacobian");
- octave_idx_type vradaujac = 0; // We need to set this if no Jac available
+ F77_INT vradaujac = 0; // We need to set this if no Jac available
if (!vradaujacfun.is_empty ()) vradaujac = 1;
// The option JPattern will be ignored by this solver, the core
@@ -493,7 +493,7 @@ ode2r (@@odepkg_equations_lorenz, [0, 25
// Implementation of the option 'Mass' has been finished, these
// options can be set by the user to another value than default
vradaumass = vodeopt.contents ("Mass");
- octave_idx_type vradaumas = 0;
+ F77_INT vradaumas = 0;
if (!vradaumass.is_empty ()) {
vradaumas = 1;
if (vradaumass.is_function_handle () || vradaumass.is_inline_function ())
@@ -567,31 +567,31 @@ ode2r (@@odepkg_equations_lorenz, [0, 25
NDArray vRTOL = vreltol.array_value ();
NDArray vATOL = vabstol.array_value ();
- octave_idx_type N = args(2).length ();
+ F77_INT N = TO_F77_INT (args(2).length ());
double X = args(1).vector_value ()(0);
double* Y = vY0.fortran_vec ();
double XEND = args(1).vector_value ()(1);
double H = vinitstep.double_value ();
double *RTOL = vRTOL.fortran_vec ();
double *ATOL = vATOL.fortran_vec ();
- octave_idx_type ITOL = vitol;
- octave_idx_type IJAC = vradaujac;
- octave_idx_type MLJAC=N;
- octave_idx_type MUJAC=N;
+ F77_INT ITOL = vitol;
+ F77_INT IJAC = vradaujac;
+ F77_INT MLJAC=N;
+ F77_INT MUJAC=N;
- octave_idx_type IMAS=vradaumas;
- octave_idx_type MLMAS=N;
- octave_idx_type MUMAS=N;
- octave_idx_type IOUT = 1; // The SOLOUT function will always be called
- octave_idx_type LWORK = N*(N+N+7*N+3*7+3)+20;
+ F77_INT IMAS=vradaumas;
+ F77_INT MLMAS=N;
+ F77_INT MUMAS=N;
+ F77_INT IOUT = 1; // The SOLOUT function will always be called
+ F77_INT LWORK = N*(N+N+7*N+3*7+3)+20;
OCTAVE_LOCAL_BUFFER (double, WORK, LWORK);
- for (octave_idx_type vcnt = 0; vcnt < LWORK; vcnt++) WORK[vcnt] = 0.0;
- octave_idx_type LIWORK = (2+(7-1)/2)*N+20;
- OCTAVE_LOCAL_BUFFER (octave_idx_type, IWORK, LIWORK);
- for (octave_idx_type vcnt = 0; vcnt < LIWORK; vcnt++) IWORK[vcnt] = 0;
+ for (F77_INT vcnt = 0; vcnt < LWORK; vcnt++) WORK[vcnt] = 0.0;
+ F77_INT LIWORK = (2+(7-1)/2)*N+20;
+ OCTAVE_LOCAL_BUFFER (F77_INT, IWORK, LIWORK);
+ for (F77_INT vcnt = 0; vcnt < LIWORK; vcnt++) IWORK[vcnt] = 0;
double RPAR[1] = {0.0};
- octave_idx_type IPAR[1] = {0};
- octave_idx_type IDID = 0;
+ F77_INT IPAR[1] = {0};
+ F77_INT IDID = 0;
IWORK[0] = 1; // Switch for transformation of Jacobian into Hessenberg form
WORK[2] = -1; // Recompute Jacobian after every succesful step
diff --git a/src/odepkg_octsolver_radau5.cc b/src/odepkg_octsolver_radau5.cc
--- a/src/odepkg_octsolver_radau5.cc
+++ b/src/odepkg_octsolver_radau5.cc
@@ -38,46 +38,46 @@ similiar to the implementation of this f
#include
#include "odepkg_auxiliary_functions.h"
-typedef octave_idx_type (*odepkg_radau5_usrtype)
- (const octave_idx_type& N, const double& X, const double* Y, double* F,
+typedef F77_RET_T (*odepkg_radau5_usrtype)
+ (const F77_INT& N, const double& X, const double* Y, double* F,
OCTAVE_UNUSED const double* RPAR,
- OCTAVE_UNUSED const octave_idx_type* IPAR);
+ OCTAVE_UNUSED const F77_INT* IPAR);
-typedef octave_idx_type (*odepkg_radau5_jactype)
- (const octave_idx_type& N, const double& X, const double* Y, double* DFY,
- OCTAVE_UNUSED const octave_idx_type& LDFY,
+typedef F77_RET_T (*odepkg_radau5_jactype)
+ (const F77_INT& N, const double& X, const double* Y, double* DFY,
+ OCTAVE_UNUSED const F77_INT& LDFY,
OCTAVE_UNUSED const double* RPAR,
- OCTAVE_UNUSED const octave_idx_type* IPAR);
+ OCTAVE_UNUSED const F77_INT* IPAR);
typedef F77_RET_T (*odepkg_radau5_masstype)
- (const octave_idx_type& N, double* AM,
- OCTAVE_UNUSED const octave_idx_type* LMAS,
+ (const F77_INT& N, double* AM,
+ OCTAVE_UNUSED const F77_INT* LMAS,
OCTAVE_UNUSED const double* RPAR,
- OCTAVE_UNUSED const octave_idx_type* IPAR);
+ OCTAVE_UNUSED const F77_INT* IPAR);
-typedef octave_idx_type (*odepkg_radau5_soltype)
- (const octave_idx_type& NR, const double& XOLD, const double& X,
- const double* Y, const double* CONT, const octave_idx_type* LRC,
- const octave_idx_type& N, OCTAVE_UNUSED const double* RPAR,
- OCTAVE_UNUSED const octave_idx_type* IPAR, octave_idx_type& IRTRN);
+typedef F77_RET_T (*odepkg_radau5_soltype)
+ (const F77_INT& NR, const double& XOLD, const double& X,
+ const double* Y, const double* CONT, const F77_INT* LRC,
+ const F77_INT& N, OCTAVE_UNUSED const double* RPAR,
+ OCTAVE_UNUSED const F77_INT* IPAR, F77_INT& IRTRN);
extern "C" {
F77_RET_T F77_FUNC (radau5, RADAU5)
- (const octave_idx_type& N, odepkg_radau5_usrtype,
+ (const F77_INT& N, odepkg_radau5_usrtype,
const double& X, const double* Y, const double& XEND,
const double& H, const double* RTOL, const double* ATOL,
- const octave_idx_type& ITOL, odepkg_radau5_jactype, const octave_idx_type& IJAC,
- const octave_idx_type& MLJAC, const octave_idx_type& MUJAC, odepkg_radau5_masstype,
- const octave_idx_type& IMAS, const octave_idx_type& MLMAS,
- const octave_idx_type& MUMAS, odepkg_radau5_soltype, const octave_idx_type& IOUT,
- const double* WORK, const octave_idx_type& LWORK, const octave_idx_type* IWORK,
- const octave_idx_type& LIWORK, OCTAVE_UNUSED const double* RPAR,
- OCTAVE_UNUSED const octave_idx_type* IPAR, const octave_idx_type& IDID);
+ const F77_INT& ITOL, odepkg_radau5_jactype, const F77_INT& IJAC,
+ const F77_INT& MLJAC, const F77_INT& MUJAC, odepkg_radau5_masstype,
+ const F77_INT& IMAS, const F77_INT& MLMAS,
+ const F77_INT& MUMAS, odepkg_radau5_soltype, const F77_INT& IOUT,
+ const double* WORK, const F77_INT& LWORK, const F77_INT* IWORK,
+ const F77_INT& LIWORK, OCTAVE_UNUSED const double* RPAR,
+ OCTAVE_UNUSED const F77_INT* IPAR, const F77_INT& IDID);
double F77_FUNC (contr5, CONTR5)
- (const octave_idx_type& I, const double& S,
- const double* CONT, const octave_idx_type* LRC);
+ (const F77_INT& I, const double& S,
+ const double* CONT, const F77_INT* LRC);
} // extern "C"
@@ -96,15 +96,15 @@ static Matrix denseOutputMatrix;
static ColumnVector denseTimeVector;
static int denseTimeIndex;
-octave_idx_type odepkg_radau5_usrfcn
- (const octave_idx_type& N, const double& X, const double* Y,
+F77_RET_T odepkg_radau5_usrfcn
+ (const F77_INT& N, const double& X, const double* Y,
double* F, OCTAVE_UNUSED const double* RPAR,
- OCTAVE_UNUSED const octave_idx_type* IPAR) {
+ OCTAVE_UNUSED const F77_INT* IPAR) {
// Copy the values that come from the Fortran function element wise,
// otherwise Octave will crash if these variables will be freed
ColumnVector A(N);
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++) {
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++) {
A(vcnt) = Y[vcnt];
// octave_stdout << "I am here Y[" << vcnt << "] " << Y[vcnt] << std::endl;
// octave_stdout << "I am here T " << X << std::endl;
@@ -114,29 +114,29 @@ octave_idx_type odepkg_radau5_usrfcn
// function that keeps the set of differential algebraic equations
octave_value_list varin;
varin(0) = X; varin(1) = A;
- for (octave_idx_type vcnt = 0; vcnt < vradau5extarg.length (); vcnt++)
+ for (F77_INT vcnt = 0; vcnt < vradau5extarg.length (); vcnt++)
varin(vcnt+2) = vradau5extarg(vcnt);
octave_value_list vout = feval (vradau5odefun.function_value (), varin, 1);
// Return the results from the function evaluation to the Fortran
// solver, again copy them and don't just create a Fortran vector
ColumnVector vcol = vout(0).column_vector_value ();
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++)
F[vcnt] = vcol(vcnt);
- return (true);
+ F77_RETURN (true);
}
-octave_idx_type odepkg_radau5_jacfcn
- (const octave_idx_type& N, const double& X, const double* Y,
- double* DFY, OCTAVE_UNUSED const octave_idx_type& LDFY,
+F77_RET_T odepkg_radau5_jacfcn
+ (const F77_INT& N, const double& X, const double* Y,
+ double* DFY, OCTAVE_UNUSED const F77_INT& LDFY,
OCTAVE_UNUSED const double* RPAR,
- OCTAVE_UNUSED const octave_idx_type* IPAR) {
+ OCTAVE_UNUSED const F77_INT* IPAR) {
// Copy the values that come from the Fortran function element-wise,
// otherwise Octave will crash if these variables are freed
ColumnVector A(N);
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++)
A(vcnt) = Y[vcnt];
// Set the values that are needed as input arguments before calling
@@ -147,23 +147,23 @@ octave_idx_type odepkg_radau5_jacfcn
(vradau5jacfun, vt, vy, vradau5extarg);
Matrix vdfy = vout.matrix_value ();
- for (octave_idx_type vcol = 0; vcol < N; vcol++)
- for (octave_idx_type vrow = 0; vrow < N; vrow++)
+ for (F77_INT vcol = 0; vcol < N; vcol++)
+ for (F77_INT vrow = 0; vrow < N; vrow++)
DFY[vrow+vcol*N] = vdfy (vrow, vcol);
- return (true);
+ F77_RETURN (true);
}
F77_RET_T odepkg_radau5_massfcn
- (const octave_idx_type& N, double* AM,
- OCTAVE_UNUSED const octave_idx_type* LMAS,
+ (const F77_INT& N, double* AM,
+ OCTAVE_UNUSED const F77_INT* LMAS,
OCTAVE_UNUSED const double* RPAR,
- OCTAVE_UNUSED const octave_idx_type* IPAR) {
+ OCTAVE_UNUSED const F77_INT* IPAR) {
// Copy the values that come from the Fortran function element-wise,
// otherwise Octave will crash if these variables are freed
ColumnVector A(N);
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++)
A(vcnt) = 0.0;
// warning_with_id ("OdePkg:InvalidFunctionCall",
@@ -177,23 +177,23 @@ F77_RET_T odepkg_radau5_massfcn
(vradau5mass, vradau5massstate, vt, vy, vradau5extarg);
Matrix vam = vout.matrix_value ();
- for (octave_idx_type vrow = 0; vrow < N; vrow++)
- for (octave_idx_type vcol = 0; vcol < N; vcol++)
+ for (F77_INT vrow = 0; vrow < N; vrow++)
+ for (F77_INT vcol = 0; vcol < N; vcol++)
AM[vrow+vcol*N] = vam (vrow, vcol);
return ((F77_RET_T) true);
}
-octave_idx_type odepkg_radau5_solfcn
- (const octave_idx_type& NR, const double& XOLD, const double& X,
- const double* Y, const double* CONT, const octave_idx_type* LRC,
- const octave_idx_type& N, OCTAVE_UNUSED const double* RPAR,
- OCTAVE_UNUSED const octave_idx_type* IPAR, octave_idx_type& IRTRN) {
+F77_RET_T odepkg_radau5_solfcn
+ (const F77_INT& NR, const double& XOLD, const double& X,
+ const double* Y, const double* CONT, const F77_INT* LRC,
+ const F77_INT& N, OCTAVE_UNUSED const double* RPAR,
+ OCTAVE_UNUSED const F77_INT* IPAR, F77_INT& IRTRN) {
// Copy the values that come from the Fortran function element-wise,
// otherwise Octave will crash if these variables are freed
ColumnVector A(N);
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++)
A(vcnt) = Y[vcnt];
// Set the values that are needed as input arguments before calling
@@ -224,12 +224,12 @@ octave_idx_type odepkg_radau5_solfcn
if (!vradau5pltfun.is_empty ()) {
if (vradau5refine.int_value () > 0) {
ColumnVector B(N); double vtb = 0.0;
- for (octave_idx_type vcnt = 1; vcnt < vradau5refine.int_value (); vcnt++) {
+ for (F77_INT vcnt = 1; vcnt < vradau5refine.int_value (); vcnt++) {
// Calculate time stamps between XOLD and X and get the
// results at these time stamps
vtb = (X - XOLD) * vcnt / vradau5refine.int_value () + XOLD;
- for (octave_idx_type vcou = 0; vcou < N; vcou++)
+ for (F77_INT vcou = 0; vcou < N; vcou++)
B(vcou) = F77_FUNC (contr5, CONTR5) (vcou+1, vtb, CONT, LRC);
// Evaluate the 'OutputFcn' with the approximated values from
@@ -242,18 +242,18 @@ octave_idx_type odepkg_radau5_solfcn
}
// Evaluate the 'OutputFcn' with the results from the solver, if
// the OutputFcn returns true then set a negative value in IRTRN
- IRTRN = - odepkg_auxiliary_evalplotfun
- (vradau5pltfun, vradau5outsel, vt, vy, vradau5extarg, 1);
+ IRTRN = TO_F77_INT (- odepkg_auxiliary_evalplotfun
+ (vradau5pltfun, vradau5outsel, vt, vy, vradau5extarg, 1));
}
// Check if dense output is required. If yes, calculate the solution at the given points
if (denseOutputStatus == true) {
- for (octave_idx_type idxDen = denseTimeIndex; idxDen < denseTimeVector.numel(); idxDen++)
+ for (F77_INT idxDen = denseTimeIndex; idxDen < denseTimeVector.numel(); idxDen++)
{
// TODO: Maybe there is a more efficient way to evaluate this
if (denseTimeVector(idxDen) >= XOLD && denseTimeVector(idxDen) <= X)
{
- for (octave_idx_type vcou = 0; vcou < N; vcou++)
+ for (F77_INT vcou = 0; vcou < N; vcou++)
{
// If value is found then evaluate with the contr5 function
denseOutputMatrix(idxDen, vcou)=F77_FUNC (contr5, CONTR5) (vcou+1, denseTimeVector(idxDen), CONT, LRC);
@@ -271,7 +271,7 @@ octave_idx_type odepkg_radau5_solfcn
// Dense output not activated ...
}
- return (true);
+ F77_RETURN (true);
}
// PKG_ADD: autoload ("ode5r", "dldsolver.oct");
@@ -427,7 +427,7 @@ ode5r (@@odepkg_equations_lorenz, [0, 25
// Setting the tolerance type that depends on the types (scalar or
// vector) of the options RelTol and AbsTol
- octave_idx_type vitol = 0;
+ F77_INT vitol = 0;
if (vreltol.is_scalar_type () && (vreltol.length () == vabstol.length ()))
vitol = 0;
else if (!vreltol.is_scalar_type () && (vreltol.length () == vabstol.length ()))
@@ -519,7 +519,7 @@ ode5r (@@odepkg_equations_lorenz, [0, 25
// Implementation of the option 'Jacobian' has been finished, these
// options can be set by the user to another value than default
vradau5jacfun = vodeopt.contents ("Jacobian");
- octave_idx_type vradau5jac = 0; // We need to set this if no Jac available
+ F77_INT vradau5jac = 0; // We need to set this if no Jac available
if (!vradau5jacfun.is_empty ()) vradau5jac = 1;
// The option JPattern will be ignored by this solver, the core
@@ -540,7 +540,7 @@ ode5r (@@odepkg_equations_lorenz, [0, 25
// Implementation of the option 'Mass' has been finished, these
// options can be set by the user to another value than default
vradau5mass = vodeopt.contents ("Mass");
- octave_idx_type vradau5mas = 0;
+ F77_INT vradau5mas = 0;
if (!vradau5mass.is_empty ()) {
vradau5mas = 1;
if (vradau5mass.is_function_handle () || vradau5mass.is_inline_function ())
@@ -623,31 +623,31 @@ ode5r (@@odepkg_equations_lorenz, [0, 25
NDArray vRTOL = vreltol.array_value ();
NDArray vATOL = vabstol.array_value ();
- octave_idx_type N = args(2).length ();
+ F77_INT N = TO_F77_INT (args(2).length ());
double X = tSpan(0);
double* Y = vY0.fortran_vec ();
double XEND = tSpan(1);
double H = vinitstep.double_value ();
double *RTOL = vRTOL.fortran_vec ();
double *ATOL = vATOL.fortran_vec ();
- octave_idx_type ITOL = vitol;
- octave_idx_type IJAC = vradau5jac;
- octave_idx_type MLJAC=N;
- octave_idx_type MUJAC=N;
+ F77_INT ITOL = vitol;
+ F77_INT IJAC = vradau5jac;
+ F77_INT MLJAC=N;
+ F77_INT MUJAC=N;
- octave_idx_type IMAS=vradau5mas;
- octave_idx_type MLMAS=N;
- octave_idx_type MUMAS=N;
- octave_idx_type IOUT = 1; // The SOLOUT function will always be called
- octave_idx_type LWORK = N*(N+N+7*N+3*7+3)+20;
+ F77_INT IMAS=vradau5mas;
+ F77_INT MLMAS=N;
+ F77_INT MUMAS=N;
+ F77_INT IOUT = 1; // The SOLOUT function will always be called
+ F77_INT LWORK = N*(N+N+7*N+3*7+3)+20;
OCTAVE_LOCAL_BUFFER (double, WORK, LWORK);
- for (octave_idx_type vcnt = 0; vcnt < LWORK; vcnt++) WORK[vcnt] = 0.0;
- octave_idx_type LIWORK = (2+(7-1)/2)*N+20;
- OCTAVE_LOCAL_BUFFER (octave_idx_type, IWORK, LIWORK);
- for (octave_idx_type vcnt = 0; vcnt < LIWORK; vcnt++) IWORK[vcnt] = 0;
+ for (F77_INT vcnt = 0; vcnt < LWORK; vcnt++) WORK[vcnt] = 0.0;
+ F77_INT LIWORK = (2+(7-1)/2)*N+20;
+ OCTAVE_LOCAL_BUFFER (F77_INT, IWORK, LIWORK);
+ for (F77_INT vcnt = 0; vcnt < LIWORK; vcnt++) IWORK[vcnt] = 0;
double RPAR[1] = {0.0};
- octave_idx_type IPAR[1] = {0};
- octave_idx_type IDID = 0;
+ F77_INT IPAR[1] = {0};
+ F77_INT IDID = 0;
IWORK[0] = 1; // Switch for transformation of Jacobian into Hessenberg form
IWORK[2] = vmaxnit.int_value (); // Maximum number of Newton iterations
diff --git a/src/odepkg_octsolver_rodas.cc b/src/odepkg_octsolver_rodas.cc
--- a/src/odepkg_octsolver_rodas.cc
+++ b/src/odepkg_octsolver_rodas.cc
@@ -36,41 +36,41 @@ Compile this file manually and run some
/* -*- texinfo -*-
* @subsection Source File @file{odepkg_octsolver_rodas.cc}
*
- * @deftp {Typedef} {octave_idx_type (*odepkg_rodas_usrtype)}
+ * @deftp {Typedef} {F77_INT (*odepkg_rodas_usrtype)}
* This @code{typedef} is used to define the input and output arguments of the user function for the DAE problem that is further needed by the Fortran core solver @code{rodas}. The implementation of this @code{typedef} is
*
* @example
- * typedef octave_idx_type (*odepkg_rodas_usrtype)
- * (const octave_idx_type& N, const double& X, const double* Y,
+ * typedef F77_RET_T (*odepkg_rodas_usrtype)
+ * (const F77_INT& N, const double& X, const double* Y,
* double* F, const double* RPAR,
- * const octave_idx_type* IPAR);
+ * const F77_INT* IPAR);
* @end example
* @end deftp
*/
-typedef octave_idx_type (*odepkg_rodas_usrtype)
- (const octave_idx_type& N, const double& X, const double* Y, double* F,
+typedef F77_RET_T (*odepkg_rodas_usrtype)
+ (const F77_INT& N, const double& X, const double* Y, double* F,
const double* RPAR,
- const octave_idx_type* IPAR);
+ const F77_INT* IPAR);
/* -*- texinfo -*-
- * @deftp {Typedef} {octave_idx_type (*odepkg_rodas_jactype)}
+ * @deftp {Typedef} {F77_INT (*odepkg_rodas_jactype)}
*
* This @code{typedef} is used to define the input and output arguments of the @code{Jacobian} function for the DAE problem that is further needed by the Fortran core solver @code{rodas}. The implementation of this @code{typedef} is
*
* @example
- * typedef octave_idx_type (*odepkg_rodas_jactype)
- * (const octave_idx_type& N, const double& X, const double* Y,
- * double* DFY, const octave_idx_type* LDFY,
+ * typedef F77_RET_T (*odepkg_rodas_jactype)
+ * (const F77_INT& N, const double& X, const double* Y,
+ * double* DFY, const F77_INT* LDFY,
* const double* RPAR,
- * const octave_idx_type* IPAR);
+ * const F77_INT* IPAR);
* @end example
* @end deftp
*/
-typedef octave_idx_type (*odepkg_rodas_jactype)
- (const octave_idx_type& N, const double& X, const double* Y, double* DFY,
- const octave_idx_type& LDFY,
+typedef F77_RET_T (*odepkg_rodas_jactype)
+ (const F77_INT& N, const double& X, const double* Y, double* DFY,
+ const F77_INT& LDFY,
const double* RPAR,
- const octave_idx_type* IPAR);
+ const F77_INT* IPAR);
/* -*- texinfo -*-
* @deftp {Typedef} {F77_RET_T (*odepkg_rodas_masstype)}
@@ -79,100 +79,100 @@ typedef octave_idx_type (*odepkg_rodas_j
*
* @example
* typedef F77_RET_T (*odepkg_rodas_masstype)
- * (const octave_idx_type& N, double* AM,
- * const octave_idx_type* LMAS,
+ * (const F77_INT& N, double* AM,
+ * const F77_INT* LMAS,
* const double* RPAR,
- * const octave_idx_type* IPAR);
+ * const F77_INT* IPAR);
* @end example
* @end deftp
*/
typedef F77_RET_T (*odepkg_rodas_masstype)
- (const octave_idx_type& N, double* AM,
- const octave_idx_type* LMAS,
+ (const F77_INT& N, double* AM,
+ const F77_INT* LMAS,
const double* RPAR,
- const octave_idx_type* IPAR);
+ const F77_INT* IPAR);
/* -*- texinfo -*-
- * @deftp {Typedef} {octave_idx_type (*odepkg_rodas_soltype)}
+ * @deftp {Typedef} {F77_INT (*odepkg_rodas_soltype)}
*
* This @code{typedef} is used to define the input and output arguments of the @code{Solution} function for the DAE problem that is further needed by the Fortran core solver @code{rodas}. The implementation of this @code{typedef} is
*
* @example
- * typedef octave_idx_type (*odepkg_rodas_soltype)
- * (const octave_idx_type& NR, const double& XOLD, const double& X,
- * const double* Y, const double* CONT, const octave_idx_type* LRC,
- * const octave_idx_type& N, const double* RPAR,
- * const octave_idx_type* IPAR, octave_idx_type& IRTRN);
+ * typedef F77_RET_T (*odepkg_rodas_soltype)
+ * (const F77_INT& NR, const double& XOLD, const double& X,
+ * const double* Y, const double* CONT, const F77_INT* LRC,
+ * const F77_INT& N, const double* RPAR,
+ * const F77_INT* IPAR, F77_INT& IRTRN);
* @end example
* @end deftp
*/
-typedef octave_idx_type (*odepkg_rodas_soltype)
- (const octave_idx_type& NR, const double& XOLD, const double& X,
- const double* Y, const double* CONT, const octave_idx_type* LRC,
- const octave_idx_type& N, const double* RPAR,
- const octave_idx_type* IPAR, octave_idx_type& IRTRN);
+typedef F77_RET_T (*odepkg_rodas_soltype)
+ (const F77_INT& NR, const double& XOLD, const double& X,
+ const double* Y, const double* CONT, const F77_INT* LRC,
+ const F77_INT& N, const double* RPAR,
+ const F77_INT* IPAR, F77_INT& IRTRN);
/* -*- texinfo -*-
- * @deftp {Typedef} {octave_idx_type (*odepkg_rodas_dfxtype)}
+ * @deftp {Typedef} {F77_INT (*odepkg_rodas_dfxtype)}
*
* This @code{typedef} is used to define the input and output arguments of the @code{DFX} function for the DAE problem that is further needed by the Fortran core solver @code{rodas}. The implementation of this @code{typedef} is
*
* @example
- * typedef octave_idx_type (*odepkg_rodas_dfxtype)
- * ( const octave_idx_type& N,
+ * typedef F77_RET_T (*odepkg_rodas_dfxtype)
+ * ( const F77_INT& N,
* const double& X,
* const double* Y,
* const double* FX,
* const double* RPAR,
- * const octave_idx_type* IPAR);
+ * const F77_INT* IPAR);
* @end example
* @end deftp
*/
-typedef octave_idx_type (*odepkg_rodas_dfxtype)
- ( const octave_idx_type& N, const double& X,
+typedef F77_RET_T (*odepkg_rodas_dfxtype)
+ ( const F77_INT& N, const double& X,
const double* Y, const double* FX,
- const double* RPAR, const octave_idx_type* IPAR);
+ const double* RPAR, const F77_INT* IPAR);
extern "C" {
/* -*- texinfo -*-
- * @deftp {Prototype} {F77_RET_T F77_FUNC (rodas, RODAS)} (const octave_idx_type& N, odepkg_rodas_usrtype, const octave_idx_type& IFCN, const octave_idx_type& X, const double* Y, const double& XEND, const double& H, const double* RTOL, const double* ATOL, const octave_idx_type& ITOL, odepkg_rodas_jactype, const octave_idx_type& IJAC, const octave_idx_type& MLJAC, const octave_idx_type& MUJAC, odepkg_rodas_dfxtype, const octave_idx_type& IDFX, odepkg_rodas_masstype, const octave_idx_type& IMAS, const octave_idx_type& MLMAS, const octave_idx_type& MUMAS, odepkg_rodas_soltype, const octave_idx_type& IOUT, const double* WORK, const octave_idx_type& LWORK, const octave_idx_type* IWORK, const octave_idx_type& LIWORK, const double* RPAR, const octave_idx_type* IPAR, const octave_idx_type& IDID);
+ * @deftp {Prototype} {F77_RET_T F77_FUNC (rodas, RODAS)} (const F77_INT& N, odepkg_rodas_usrtype, const F77_INT& IFCN, const F77_INT& X, const double* Y, const double& XEND, const double& H, const double* RTOL, const double* ATOL, const F77_INT& ITOL, odepkg_rodas_jactype, const F77_INT& IJAC, const F77_INT& MLJAC, const F77_INT& MUJAC, odepkg_rodas_dfxtype, const F77_INT& IDFX, odepkg_rodas_masstype, const F77_INT& IMAS, const F77_INT& MLMAS, const F77_INT& MUMAS, odepkg_rodas_soltype, const F77_INT& IOUT, const double* WORK, const F77_INT& LWORK, const F77_INT* IWORK, const F77_INT& LIWORK, const double* RPAR, const F77_INT* IPAR, const F77_INT& IDID);
*
* The prototype @code{F77_FUNC (rodas, RODAS)} is used to represent the information about the Fortran core solver @code{rodas} that is defined in the Fortran source file @file{rodas.f} (cf. the Fortran source file @file{rodas.f} for further details).
* @end deftp
*/
F77_RET_T F77_FUNC (rodas, RODAS)
- (const octave_idx_type& N, odepkg_rodas_usrtype, const octave_idx_type& IFCN,
+ (const F77_INT& N, odepkg_rodas_usrtype, const F77_INT& IFCN,
const double& X, const double* Y, const double& XEND,
const double& H, const double* RTOL, const double* ATOL,
- const octave_idx_type& ITOL, odepkg_rodas_jactype, const octave_idx_type& IJAC,
- const octave_idx_type& MLJAC, const octave_idx_type& MUJAC, odepkg_rodas_dfxtype,
- const octave_idx_type& IDFX, odepkg_rodas_masstype, const octave_idx_type& IMAS,
- const octave_idx_type& MLMAS, const octave_idx_type& MUMAS, odepkg_rodas_soltype,
- const octave_idx_type& IOUT, const double* WORK, const octave_idx_type& LWORK,
- const octave_idx_type* IWORK, const octave_idx_type& LIWORK,
- const double* RPAR, const octave_idx_type* IPAR,
- const octave_idx_type& IDID);
+ const F77_INT& ITOL, odepkg_rodas_jactype, const F77_INT& IJAC,
+ const F77_INT& MLJAC, const F77_INT& MUJAC, odepkg_rodas_dfxtype,
+ const F77_INT& IDFX, odepkg_rodas_masstype, const F77_INT& IMAS,
+ const F77_INT& MLMAS, const F77_INT& MUMAS, odepkg_rodas_soltype,
+ const F77_INT& IOUT, const double* WORK, const F77_INT& LWORK,
+ const F77_INT* IWORK, const F77_INT& LIWORK,
+ const double* RPAR, const F77_INT* IPAR,
+ const F77_INT& IDID);
/* -*- texinfo -*-
- * @deftp {Prototype} {double F77_FUNC (contro, CONTRO)} (const octave_idx_type& I, const double* S, const double* CONT, const octave_idx_type* LRC);
+ * @deftp {Prototype} {double F77_FUNC (contro, CONTRO)} (const F77_INT& I, const double* S, const double* CONT, const F77_INT* LRC);
*
* The prototype @code{F77_FUNC (contro, CONTRO)} is used to represent the information about a continous output calculation for the @code{rodas} algorithm that is defined in the Fortran source file @file{rodas.f} (cf. the Fortran source file @file{rodas.f} for further details).
* @end deftp
*/
double F77_FUNC (contro, CONTRO)
- (const octave_idx_type& I, const double& S,
- const double* CONT, const octave_idx_type* LRC);
+ (const F77_INT& I, const double& S,
+ const double* CONT, const F77_INT* LRC);
/* -*- texinfo -*-
- * @deftp {Prototype} {F77_RET_T F77_FUNC (dfx, DFX)} ( const octave_idx_type& N, const double& X, const double* Y, const double* FX, const double* RPAR, const octave_idx_type* IPAR);
+ * @deftp {Prototype} {F77_RET_T F77_FUNC (dfx, DFX)} ( const F77_INT& N, const double& X, const double* Y, const double* FX, const double* RPAR, const F77_INT* IPAR);
*
* The prototype @code{F77_RET_T F77_FUNC (dfx, DFX)} is used to represent the information about a DFX dummy function (that will never be called because it is not supported) for the @code{rodas} algorithm that is defined in the Fortran source file @file{rodas.f} (cf. the Fortran source file @file{rodas.f} for further details).
* @end deftp
*/
F77_RET_T F77_FUNC (dfx, DFX)
- ( const octave_idx_type& N, const double& X,
+ ( const F77_INT& N, const double& X,
const double* Y, const double* FX,
- const double* RPAR, const octave_idx_type* IPAR);
+ const double* RPAR, const F77_INT* IPAR);
} // extern "C"
/* -*- texinfo -*-
@@ -256,7 +256,7 @@ static octave_value vrodasmass;
static octave_value vrodasmassstate;
/* -*- texinfo -*-
- * @deftypefn {Function} {octave_idx_type} {odepkg_rodas_usrfcn} (const octave_idx_type& N, const double& X, const double* Y, double* F, const double* RPAR, const octave_idx_type* IPAR)
+ * @deftypefn {Function} {F77_INT} {odepkg_rodas_usrfcn} (const F77_INT& N, const double& X, const double* Y, double* F, const double* RPAR, const F77_INT* IPAR)
*
* Return @code{true} if the evaluation of the user type function (that is defined for a special DAE problem in Octave) was successful, otherwise return @code{false}. This function is directly called from the Fortran core solver @code{rodas}. The input arguments of this function are
*
@@ -270,15 +270,15 @@ static octave_value vrodasmassstate;
* @end itemize
* @end deftypefn
*/
-octave_idx_type odepkg_rodas_usrfcn
- (const octave_idx_type& N, const double& X, const double* Y,
+F77_RET_T odepkg_rodas_usrfcn
+ (const F77_INT& N, const double& X, const double* Y,
double* F, const double* RPAR,
- const octave_idx_type* IPAR) {
+ const F77_INT* IPAR) {
// Copy the values that come from the Fortran function element wise,
// otherwise Octave will crash if these variables will be freed
ColumnVector A(N);
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++) {
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++) {
A(vcnt) = Y[vcnt];
// octave_stdout << "I am here Y[" << vcnt << "] " << Y[vcnt] << std::endl;
// octave_stdout << "I am here T " << X << std::endl;
@@ -288,21 +288,21 @@ octave_idx_type odepkg_rodas_usrfcn
// function that keeps the set of differential algebraic equations
octave_value_list varin;
varin(0) = X; varin(1) = A;
- for (octave_idx_type vcnt = 0; vcnt < vrodasextarg.length (); vcnt++)
+ for (F77_INT vcnt = 0; vcnt < vrodasextarg.length (); vcnt++)
varin(vcnt+2) = vrodasextarg(vcnt);
octave_value_list vout = feval (vrodasodefun.function_value (), varin, 1);
// Return the results from the function evaluation to the Fortran
// solver, again copy them and don't just create a Fortran vector
ColumnVector vcol = vout(0).column_vector_value ();
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++)
F[vcnt] = vcol(vcnt);
- return (true);
+ F77_RETURN (true);
}
/* -*- texinfo -*-
- * @deftypefn {Function} {octave_idx_type} {odepkg_mebdfi_jacfcn} (const octave_idx_type& N, const double& X, const double* Y, double* DFY, const octave_idx_type& LDFY, const double* RPAR, const octave_idx_type* IPAR)
+ * @deftypefn {Function} {F77_INT} {odepkg_mebdfi_jacfcn} (const F77_INT& N, const double& X, const double* Y, double* DFY, const F77_INT& LDFY, const double* RPAR, const F77_INT* IPAR)
*
* Return @code{true} if the evaluation of the @code{Jacobian} function (that is defined for a special DAE problem in Octave) was successful, otherwise return @code{false}. This function is directly called from the Fortran core solver @code{rodas}. The input arguments of this function are
*
@@ -317,16 +317,16 @@ octave_idx_type odepkg_rodas_usrfcn
* @end itemize
* @end deftypefn
*/
-octave_idx_type odepkg_rodas_jacfcn
- (const octave_idx_type& N, const double& X, const double* Y,
- double* DFY, const octave_idx_type& LDFY,
+F77_RET_T odepkg_rodas_jacfcn
+ (const F77_INT& N, const double& X, const double* Y,
+ double* DFY, const F77_INT& LDFY,
const double* RPAR,
- const octave_idx_type* IPAR) {
+ const F77_INT* IPAR) {
// Copy the values that come from the Fortran function element-wise,
// otherwise Octave will crash if these variables are freed
ColumnVector A(N);
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++)
A(vcnt) = Y[vcnt];
// Set the values that are needed as input arguments before calling
@@ -337,26 +337,26 @@ octave_idx_type odepkg_rodas_jacfcn
(vrodasjacfun, vt, vy, vrodasextarg);
Matrix vdfy = vout.matrix_value ();
- for (octave_idx_type vcol = 0; vcol < N; vcol++)
- for (octave_idx_type vrow = 0; vrow < N; vrow++)
+ for (F77_INT vcol = 0; vcol < N; vcol++)
+ for (F77_INT vrow = 0; vrow < N; vrow++)
DFY[vrow+vcol*N] = vdfy (vrow, vcol);
- return (true);
+ F77_RETURN (true);
}
/* -*- texinfo -*-
* odepkg_rodas_massfcn - TODO
*/
F77_RET_T odepkg_rodas_massfcn
- (const octave_idx_type& N, double* AM,
- const octave_idx_type* LMAS,
+ (const F77_INT& N, double* AM,
+ const F77_INT* LMAS,
const double* RPAR,
- const octave_idx_type* IPAR) {
+ const F77_INT* IPAR) {
// Copy the values that come from the Fortran function element-wise,
// otherwise Octave will crash if these variables are freed
ColumnVector A(N);
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++)
A(vcnt) = 0.0;
// warning_with_id ("OdePkg:InvalidFunctionCall",
@@ -370,8 +370,8 @@ F77_RET_T odepkg_rodas_massfcn
(vrodasmass, vrodasmassstate, vt, vy, vrodasextarg);
Matrix vam = vout.matrix_value ();
- for (octave_idx_type vrow = 0; vrow < N; vrow++)
- for (octave_idx_type vcol = 0; vcol < N; vcol++)
+ for (F77_INT vrow = 0; vrow < N; vrow++)
+ for (F77_INT vcol = 0; vcol < N; vcol++)
AM[vrow+vcol*N] = vam (vrow, vcol);
return ((F77_RET_T) true);
@@ -380,16 +380,16 @@ F77_RET_T odepkg_rodas_massfcn
/* -*- texinfo -*-
* odepkg_rodas_solfcn - TODO
*/
-octave_idx_type odepkg_rodas_solfcn
- (const octave_idx_type& NR, const double& XOLD, const double& X,
- const double* Y, const double* CONT, const octave_idx_type* LRC,
- const octave_idx_type& N, const double* RPAR,
- const octave_idx_type* IPAR, octave_idx_type& IRTRN) {
+F77_RET_T odepkg_rodas_solfcn
+ (const F77_INT& NR, const double& XOLD, const double& X,
+ const double* Y, const double* CONT, const F77_INT* LRC,
+ const F77_INT& N, const double* RPAR,
+ const F77_INT* IPAR, F77_INT& IRTRN) {
// Copy the values that come from the Fortran function element-wise,
// otherwise Octave will crash if these variables are freed
ColumnVector A(N);
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++)
A(vcnt) = Y[vcnt];
// Set the values that are needed as input arguments before calling
@@ -420,12 +420,12 @@ octave_idx_type odepkg_rodas_solfcn
if (!vrodaspltfun.is_empty ()) {
if (vrodasrefine.int_value () > 0) {
ColumnVector B(N); double vtb = 0.0;
- for (octave_idx_type vcnt = 1; vcnt < vrodasrefine.int_value (); vcnt++) {
+ for (F77_INT vcnt = 1; vcnt < vrodasrefine.int_value (); vcnt++) {
// Calculate time stamps between XOLD and X and get the
// results at these time stamps
vtb = (X - XOLD) * vcnt / vrodasrefine.int_value () + XOLD;
- for (octave_idx_type vcou = 0; vcou < N; vcou++)
+ for (F77_INT vcou = 0; vcou < N; vcou++)
B(vcou) = F77_FUNC (contro, CONTRO) (vcou+1, vtb, CONT, LRC);
// Evaluate the 'OutputFcn' with the approximated values from
@@ -438,24 +438,24 @@ octave_idx_type odepkg_rodas_solfcn
}
// Evaluate the 'OutputFcn' with the results from the solver, if
// the OutputFcn returns true then set a negative value in IRTRN
- IRTRN = - odepkg_auxiliary_evalplotfun
- (vrodaspltfun, vrodasoutsel, vt, vy, vrodasextarg, 1);
+ IRTRN = TO_F77_INT (- odepkg_auxiliary_evalplotfun
+ (vrodaspltfun, vrodasoutsel, vt, vy, vrodasextarg, 1));
}
- return (true);
+ F77_RETURN (true);
}
/* -*- texinfo -*-
* odepkg_rodas_solfcn - TODO dummy function
*/
-octave_idx_type odepkg_rodas_dfxfcn
- ( const octave_idx_type& N, const double& X,
+F77_RET_T odepkg_rodas_dfxfcn
+ ( const F77_INT& N, const double& X,
const double* Y, const double* FX,
- const double* RPAR, const octave_idx_type* IPAR) {
+ const double* RPAR, const F77_INT* IPAR) {
warning_with_id ("OdePkg:InvalidFunctionCall",
"function odepkg_rodas_dfxfcn: This warning message should never appear");
- return (true);
+ F77_RETURN (true);
}
// PKG_ADD: autoload ("oders", "dldsolver.oct");
@@ -611,7 +611,7 @@ oders (@@odepkg_equations_lorenz, [0, 25
// Setting the tolerance type that depends on the types (scalar or
// vector) of the options RelTol and AbsTol
- octave_idx_type vitol = 0;
+ F77_INT vitol = 0;
if (vreltol.is_scalar_type () && (vreltol.length () == vabstol.length ()))
vitol = 0;
else if (!vreltol.is_scalar_type () && (vreltol.length () == vabstol.length ()))
@@ -686,7 +686,7 @@ oders (@@odepkg_equations_lorenz, [0, 25
// Implementation of the option 'Jacobian' has been finished, these
// options can be set by the user to another value than default
vrodasjacfun = vodeopt.contents ("Jacobian");
- octave_idx_type vrodasjac = 0; // We need to set this if no Jac available
+ F77_INT vrodasjac = 0; // We need to set this if no Jac available
if (!vrodasjacfun.is_empty ()) vrodasjac = 1;
// The option JPattern will be ignored by this solver, the core
@@ -707,7 +707,7 @@ oders (@@odepkg_equations_lorenz, [0, 25
// Implementation of the option 'Mass' has been finished, these
// options can be set by the user to another value than default
vrodasmass = vodeopt.contents ("Mass");
- octave_idx_type vrodasmas = 0;
+ F77_INT vrodasmas = 0;
if (!vrodasmass.is_empty ()) {
vrodasmas = 1;
if (vrodasmass.is_function_handle () || vrodasmass.is_inline_function ())
@@ -780,33 +780,33 @@ oders (@@odepkg_equations_lorenz, [0, 25
NDArray vRTOL = vreltol.array_value ();
NDArray vATOL = vabstol.array_value ();
- octave_idx_type N = args(2).length ();
- octave_idx_type IFCN = 1;
+ F77_INT N = TO_F77_INT (args(2).length ());
+ F77_INT IFCN = 1;
double X = args(1).vector_value ()(0);
double* Y = vY0.fortran_vec ();
double XEND = args(1).vector_value ()(1);
double H = vinitstep.double_value ();
double *RTOL = vRTOL.fortran_vec ();
double *ATOL = vATOL.fortran_vec ();
- octave_idx_type ITOL = vitol;
- octave_idx_type IJAC = vrodasjac;
- octave_idx_type MLJAC=N;
- octave_idx_type MUJAC=N;
+ F77_INT ITOL = vitol;
+ F77_INT IJAC = vrodasjac;
+ F77_INT MLJAC=N;
+ F77_INT MUJAC=N;
- octave_idx_type IDFX=0;
- octave_idx_type IMAS=vrodasmas;
- octave_idx_type MLMAS=N;
- octave_idx_type MUMAS=N;
- octave_idx_type IOUT = 1; // The SOLOUT function will always be called
- octave_idx_type LWORK = N*(N+N+N+14)+20;
+ F77_INT IDFX=0;
+ F77_INT IMAS=vrodasmas;
+ F77_INT MLMAS=N;
+ F77_INT MUMAS=N;
+ F77_INT IOUT = 1; // The SOLOUT function will always be called
+ F77_INT LWORK = N*(N+N+N+14)+20;
OCTAVE_LOCAL_BUFFER (double, WORK, LWORK);
- for (octave_idx_type vcnt = 0; vcnt < LWORK; vcnt++) WORK[vcnt] = 0.0;
- octave_idx_type LIWORK = N+20;
- OCTAVE_LOCAL_BUFFER (octave_idx_type, IWORK, LIWORK);
- for (octave_idx_type vcnt = 0; vcnt < LIWORK; vcnt++) IWORK[vcnt] = 0;
+ for (F77_INT vcnt = 0; vcnt < LWORK; vcnt++) WORK[vcnt] = 0.0;
+ F77_INT LIWORK = N+20;
+ OCTAVE_LOCAL_BUFFER (F77_INT, IWORK, LIWORK);
+ for (F77_INT vcnt = 0; vcnt < LIWORK; vcnt++) IWORK[vcnt] = 0;
double RPAR[1] = {0.0};
- octave_idx_type IPAR[1] = {0};
- octave_idx_type IDID = 0;
+ F77_INT IPAR[1] = {0};
+ F77_INT IDID = 0;
WORK[1] = vmaxstep.double_value (); // Set the maximum step size
diff --git a/src/odepkg_octsolver_seulex.cc b/src/odepkg_octsolver_seulex.cc
--- a/src/odepkg_octsolver_seulex.cc
+++ b/src/odepkg_octsolver_seulex.cc
@@ -37,47 +37,47 @@ similiar to the implementation of this f
#include
#include "odepkg_auxiliary_functions.h"
-typedef octave_idx_type (*odepkg_seulex_usrtype)
- (const octave_idx_type& N, const double& X, const double* Y,
+typedef F77_RET_T (*odepkg_seulex_usrtype)
+ (const F77_INT& N, const double& X, const double* Y,
double* F, const double* RPAR,
- const octave_idx_type* IPAR);
+ const F77_INT* IPAR);
-typedef octave_idx_type (*odepkg_seulex_jactype)
- (const octave_idx_type& N, const double& X, const double* Y,
- double* DFY, const octave_idx_type& LDFY,
+typedef F77_RET_T (*odepkg_seulex_jactype)
+ (const F77_INT& N, const double& X, const double* Y,
+ double* DFY, const F77_INT& LDFY,
const double* RPAR,
- const octave_idx_type* IPAR);
+ const F77_INT* IPAR);
typedef F77_RET_T (*odepkg_seulex_masstype)
- (const octave_idx_type& N, double* AM,
- const octave_idx_type* LMAS,
+ (const F77_INT& N, double* AM,
+ const F77_INT* LMAS,
const double* RPAR,
- const octave_idx_type* IPAR);
+ const F77_INT* IPAR);
-typedef octave_idx_type (*odepkg_seulex_soltype)
- (const octave_idx_type& NR, const double& XOLD, const double& X,
- const double* Y, const double* RC, const octave_idx_type& LRC,
- const double* IC, const octave_idx_type& LIC,
- const octave_idx_type& N, const double* RPAR,
- const octave_idx_type* IPAR, octave_idx_type& IRTRN);
+typedef F77_RET_T (*odepkg_seulex_soltype)
+ (const F77_INT& NR, const double& XOLD, const double& X,
+ const double* Y, const double* RC, const F77_INT& LRC,
+ const double* IC, const F77_INT& LIC,
+ const F77_INT& N, const double* RPAR,
+ const F77_INT* IPAR, F77_INT& IRTRN);
extern "C" {
F77_RET_T F77_FUNC (seulex, SEULEX)
- (const octave_idx_type& N, odepkg_seulex_usrtype, const octave_idx_type& IFCN,
+ (const F77_INT& N, odepkg_seulex_usrtype, const F77_INT& IFCN,
const double& X, const double* Y, const double& XEND,
const double& H, const double* RTOL, const double* ATOL,
- const octave_idx_type& ITOL, odepkg_seulex_jactype, const octave_idx_type& IJAC,
- const octave_idx_type& MLJAC, const octave_idx_type& MUJAC, odepkg_seulex_masstype,
- const octave_idx_type& IMAS, const octave_idx_type& MLMAS,
- const octave_idx_type& MUMAS, odepkg_seulex_soltype, const octave_idx_type& IOUT,
- const double* WORK, const octave_idx_type& LWORK, const octave_idx_type* IWORK,
- const octave_idx_type& LIWORK, const double* RPAR,
- const octave_idx_type* IPAR, const octave_idx_type& IDID);
+ const F77_INT& ITOL, odepkg_seulex_jactype, const F77_INT& IJAC,
+ const F77_INT& MLJAC, const F77_INT& MUJAC, odepkg_seulex_masstype,
+ const F77_INT& IMAS, const F77_INT& MLMAS,
+ const F77_INT& MUMAS, odepkg_seulex_soltype, const F77_INT& IOUT,
+ const double* WORK, const F77_INT& LWORK, const F77_INT* IWORK,
+ const F77_INT& LIWORK, const double* RPAR,
+ const F77_INT* IPAR, const F77_INT& IDID);
double F77_FUNC (contex, CONTEX)
- (const octave_idx_type& I, const double& S,
- const double* RC, const octave_idx_type& LRC,
- const double* IC, const octave_idx_type& LIC);
+ (const F77_INT& I, const double& S,
+ const double* RC, const F77_INT& LRC,
+ const double* IC, const F77_INT& LIC);
} // extern "C"
@@ -94,15 +94,15 @@ static octave_value vseulexrefine;
static octave_value vseulexmass;
static octave_value vseulexmassstate;
-octave_idx_type odepkg_seulex_usrfcn
- (const octave_idx_type& N, const double& X, const double* Y,
+F77_RET_T odepkg_seulex_usrfcn
+ (const F77_INT& N, const double& X, const double* Y,
double* F, const double* RPAR,
- const octave_idx_type* IPAR) {
+ const F77_INT* IPAR) {
// Copy the values that come from the Fortran function element wise,
// otherwise Octave will crash if these variables will be freed
ColumnVector A(N);
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++) {
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++) {
A(vcnt) = Y[vcnt];
// octave_stdout << "I am here Y[" << vcnt << "] " << Y[vcnt] << std::endl;
// octave_stdout << "I am here T " << X << std::endl;
@@ -112,29 +112,29 @@ octave_idx_type odepkg_seulex_usrfcn
// function that keeps the set of differential algebraic equations
octave_value_list varin;
varin(0) = X; varin(1) = A;
- for (octave_idx_type vcnt = 0; vcnt < vseulexextarg.length (); vcnt++)
+ for (F77_INT vcnt = 0; vcnt < vseulexextarg.length (); vcnt++)
varin(vcnt+2) = vseulexextarg(vcnt);
octave_value_list vout = feval (vseulexodefun.function_value (), varin, 1);
// Return the results from the function evaluation to the Fortran
// solver, again copy them and don't just create a Fortran vector
ColumnVector vcol = vout(0).column_vector_value ();
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++)
F[vcnt] = vcol(vcnt);
- return (true);
+ F77_RETURN (true);
}
-octave_idx_type odepkg_seulex_jacfcn
- (const octave_idx_type& N, const double& X, const double* Y,
- double* DFY, const octave_idx_type& LDFY,
+F77_RET_T odepkg_seulex_jacfcn
+ (const F77_INT& N, const double& X, const double* Y,
+ double* DFY, const F77_INT& LDFY,
const double* RPAR,
- const octave_idx_type* IPAR) {
+ const F77_INT* IPAR) {
// Copy the values that come from the Fortran function element-wise,
// otherwise Octave will crash if these variables are freed
ColumnVector A(N);
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++)
A(vcnt) = Y[vcnt];
// Set the values that are needed as input arguments before calling
@@ -145,23 +145,23 @@ octave_idx_type odepkg_seulex_jacfcn
(vseulexjacfun, vt, vy, vseulexextarg);
Matrix vdfy = vout.matrix_value ();
- for (octave_idx_type vcol = 0; vcol < N; vcol++)
- for (octave_idx_type vrow = 0; vrow < N; vrow++)
+ for (F77_INT vcol = 0; vcol < N; vcol++)
+ for (F77_INT vrow = 0; vrow < N; vrow++)
DFY[vrow+vcol*N] = vdfy (vrow, vcol);
- return (true);
+ F77_RETURN (true);
}
F77_RET_T odepkg_seulex_massfcn
- (const octave_idx_type& N, double* AM,
- const octave_idx_type* LMAS,
+ (const F77_INT& N, double* AM,
+ const F77_INT* LMAS,
const double* RPAR,
- const octave_idx_type* IPAR) {
+ const F77_INT* IPAR) {
// Copy the values that come from the Fortran function element-wise,
// otherwise Octave will crash if these variables are freed
ColumnVector A(N);
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++)
A(vcnt) = 0.0;
// warning_with_id ("OdePkg:InvalidFunctionCall",
@@ -175,24 +175,24 @@ F77_RET_T odepkg_seulex_massfcn
(vseulexmass, vseulexmassstate, vt, vy, vseulexextarg);
Matrix vam = vout.matrix_value ();
- for (octave_idx_type vrow = 0; vrow < N; vrow++)
- for (octave_idx_type vcol = 0; vcol < N; vcol++)
+ for (F77_INT vrow = 0; vrow < N; vrow++)
+ for (F77_INT vcol = 0; vcol < N; vcol++)
AM[vrow+vcol*N] = vam (vrow, vcol);
return ((F77_RET_T) true);
}
-octave_idx_type odepkg_seulex_solfcn
- (const octave_idx_type& NR, const double& XOLD, const double& X,
- const double* Y, const double* RC, const octave_idx_type& LRC,
- const double* IC, const octave_idx_type& LIC,
- const octave_idx_type& N, const double* RPAR,
- const octave_idx_type* IPAR, octave_idx_type& IRTRN) {
+F77_RET_T odepkg_seulex_solfcn
+ (const F77_INT& NR, const double& XOLD, const double& X,
+ const double* Y, const double* RC, const F77_INT& LRC,
+ const double* IC, const F77_INT& LIC,
+ const F77_INT& N, const double* RPAR,
+ const F77_INT* IPAR, F77_INT& IRTRN) {
// Copy the values that come from the Fortran function element-wise,
// otherwise Octave will crash if these variables are freed
ColumnVector A(N);
- for (octave_idx_type vcnt = 0; vcnt < N; vcnt++)
+ for (F77_INT vcnt = 0; vcnt < N; vcnt++)
A(vcnt) = Y[vcnt];
// Set the values that are needed as input arguments before calling
@@ -226,12 +226,12 @@ octave_idx_type odepkg_seulex_solfcn
if (!vseulexpltfun.is_empty ()) {
if (vseulexrefine.int_value () > 0) {
ColumnVector B(N); double vtb = 0.0;
- for (octave_idx_type vcnt = 1; vcnt < vseulexrefine.int_value (); vcnt++) {
+ for (F77_INT vcnt = 1; vcnt < vseulexrefine.int_value (); vcnt++) {
// Calculate time stamps between XOLD and X and get the
// results at these time stamps
vtb = (X - XOLD) * vcnt / vseulexrefine.int_value () + XOLD;
- for (octave_idx_type vcou = 0; vcou < N; vcou++)
+ for (F77_INT vcou = 0; vcou < N; vcou++)
B(vcou) = F77_FUNC (contex, CONTEX) (vcou+1, vtb, RC, LRC, IC, LIC);
// Evaluate the 'OutputFcn' with the approximated values from
@@ -244,12 +244,12 @@ octave_idx_type odepkg_seulex_solfcn
}
// Evaluate the 'OutputFcn' with the results from the solver, if
// the OutputFcn returns true then set a negative value in IRTRN
- IRTRN = - odepkg_auxiliary_evalplotfun
- (vseulexpltfun, vseulexoutsel, vt, vy, vseulexextarg, 1);
+ IRTRN = TO_F77_INT (- odepkg_auxiliary_evalplotfun
+ (vseulexpltfun, vseulexoutsel, vt, vy, vseulexextarg, 1));
vseulexpltbrk = true;
}
- return (true);
+ F77_RETURN (true);
}
// PKG_ADD: autoload ("odesx", "dldsolver.oct");
@@ -405,7 +405,7 @@ odesx (@@odepkg_equations_lorenz, [0, 25
// Setting the tolerance type that depends on the types (scalar or
// vector) of the options RelTol and AbsTol
- octave_idx_type vitol = 0;
+ F77_INT vitol = 0;
if (vreltol.is_scalar_type () && (vreltol.length () == vabstol.length ()))
vitol = 0;
else if (!vreltol.is_scalar_type () && (vreltol.length () == vabstol.length ()))
@@ -480,7 +480,7 @@ odesx (@@odepkg_equations_lorenz, [0, 25
// Implementation of the option 'Jacobian' has been finished, these
// options can be set by the user to another value than default
vseulexjacfun = vodeopt.contents ("Jacobian");
- octave_idx_type vseulexjac = 0; // We need to set this if no Jac available
+ F77_INT vseulexjac = 0; // We need to set this if no Jac available
if (!vseulexjacfun.is_empty ()) vseulexjac = 1;
// The option JPattern will be ignored by this solver, the core
@@ -501,7 +501,7 @@ odesx (@@odepkg_equations_lorenz, [0, 25
// Implementation of the option 'Mass' has been finished, these
// options can be set by the user to another value than default
vseulexmass = vodeopt.contents ("Mass");
- octave_idx_type vseulexmas = 0;
+ F77_INT vseulexmas = 0;
if (!vseulexmass.is_empty ()) {
vseulexmas = 1;
if (vseulexmass.is_function_handle () || vseulexmass.is_inline_function ())
@@ -574,32 +574,32 @@ odesx (@@odepkg_equations_lorenz, [0, 25
NDArray vRTOL = vreltol.array_value ();
NDArray vATOL = vabstol.array_value ();
- octave_idx_type N = args(2).length ();
- octave_idx_type IFCN = 1;
+ F77_INT N = TO_F77_INT (args(2).length ());
+ F77_INT IFCN = 1;
double X = args(1).vector_value ()(0);
double* Y = vY0.fortran_vec ();
double XEND = args(1).vector_value ()(1);
double H = vinitstep.double_value ();
double *RTOL = vRTOL.fortran_vec ();
double *ATOL = vATOL.fortran_vec ();
- octave_idx_type ITOL = vitol;
- octave_idx_type IJAC = vseulexjac;
- octave_idx_type MLJAC=N;
- octave_idx_type MUJAC=N;
+ F77_INT ITOL = vitol;
+ F77_INT IJAC = vseulexjac;
+ F77_INT MLJAC=N;
+ F77_INT MUJAC=N;
- octave_idx_type IMAS=vseulexmas;
- octave_idx_type MLMAS=N;
- octave_idx_type MUMAS=N;
- octave_idx_type IOUT = 1; // The SOLOUT function will always be called
- octave_idx_type LWORK = N*(N+N+N+12+8)+4*12+20+(2+12*(12+3)/2)*N;
+ F77_INT IMAS=vseulexmas;
+ F77_INT MLMAS=N;
+ F77_INT MUMAS=N;
+ F77_INT IOUT = 1; // The SOLOUT function will always be called
+ F77_INT LWORK = N*(N+N+N+12+8)+4*12+20+(2+12*(12+3)/2)*N;
OCTAVE_LOCAL_BUFFER (double, WORK, LWORK);
- for (octave_idx_type vcnt = 0; vcnt < LWORK; vcnt++) WORK[vcnt] = 0.0;
- octave_idx_type LIWORK = 2*N+12+20+N;
- OCTAVE_LOCAL_BUFFER (octave_idx_type, IWORK, LIWORK);
- for (octave_idx_type vcnt = 0; vcnt < LIWORK; vcnt++) IWORK[vcnt] = 0;
+ for (F77_INT vcnt = 0; vcnt < LWORK; vcnt++) WORK[vcnt] = 0.0;
+ F77_INT LIWORK = 2*N+12+20+N;
+ OCTAVE_LOCAL_BUFFER (F77_INT, IWORK, LIWORK);
+ for (F77_INT vcnt = 0; vcnt < LIWORK; vcnt++) IWORK[vcnt] = 0;
double RPAR[1] = {0.0};
- octave_idx_type IPAR[1] = {0};
- octave_idx_type IDID = 0;
+ F77_INT IPAR[1] = {0};
+ F77_INT IDID = 0;
IWORK[0] = 1; // Switch for transformation of Jacobian into Hessenberg form
WORK[2] = -1; // Recompute Jacobian after every succesful step