# 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