/*
Copyright (C) 2000 Marc Compere, compere@mail.utexas.edu

All Octave library functions and compiler macros were copied from
examples found in lsode.cc.
The one-line copyright notice for lsode.cc is included below:

Copyright (C) 1996, 1997 John W. Eaton

The file sdirk4.cc is intended for use with Octave.

This program is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2, or (at your option) any
later version.

This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received a copy of the GNU General Public License
along with Octave; see the file COPYING.  If not, write to the Free
Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

*/
#if defined (__GNUG__)
#pragma implementation
#endif

#include "sdirk4.h"


int sdirk4_user_function (int* nstates, double* t,
                           double* x, double* retval)
{
  args.resize(2);
  args(1) = *t;

  ColumnVector V(*nstates);
  for (int i = 0; i < *nstates; i++)
    V(i) = x[i];   // <-- this is the transition from c++ to octave
  octave_value state(V);
  args(0) = state;

  if (sdirk4_fcn)
    {
      // this line is what evaluates the user-defined function
      // the stuff before and after are all lines to convert
      // fout and fargs into comaptible argument types for the
      //  cost_fcn->do_index_op command.
      // both fout and fargs are octave_value_list's
      fout = sdirk4_fcn->do_index_op( 1, args ) ;

      if (error_state)
	{
	  gripe_user_supplied_eval ("sdirk4: sdirk4_fcn, error_state");
	  return 0;
	}

      if (fout.length () > 0 && fout(0).is_defined ())
      {
	ColumnVector tmp(fout(0).vector_value()) ;
	for (int i=0; i < tmp.capacity() ; i++)
        {
          retval[i]=tmp(i);
          //fprintf( stderr, "sdirk4.cc: x[%d]=%f\n",i,tmp(i));
        }
      }

        if (error_state || fout.length () == 0)
          gripe_user_supplied_eval ("sdirk4: sdirk4_fcn, is_defined");
    }

  return 1;
}

#define SDIRK4_ABORT() \
  do \
    { \
      unwind_protect::run_frame ("Fsdirk4"); \
      return retval; \
    } \
  while (0)
 
#define SDIRK4_ABORT1(msg) \
  do \
    { \
      ::error ("sdirk4: " ## msg); \
      SDIRK4_ABORT (); \
    } \
  while (0)

#define SDIRK4_ABORT2(fmt, arg) \
  do \
    { \
      ::error ("sdirk4: " ## fmt, arg); \
      SDIRK4_ABORT (); \
    } \
  while (0)

DEFUN_DLD (sdirk4, args, nargout,
  "-*- texinfo -*-\n\
@deftypefn {Loadable Function} [@var{xout},@var{yout}] = {} sdirk4 (@var{fcn},\n\
@var{y0}, @var{xspan}, <@var{rel_tol}>, <@var{abs_tol}>, <@var{h_initial}>,\n\
<@var{abscissa_trace}>, <@var{skip_step}>)\n\
\n\
sdirk4() returns matricies of @var{xout} and @var{yout} as the solution to the\n\
ordinary differential equations\n\
specified by @var{fcn} whose initial conditions are @var{y0}.\n\
Each row in the output matrix @var{yout} corresponds\n\
to one of the elements in the vector @var{xout}.  The first element of\n\
@var{x} corresponds to the initial state @var{y0}, so that the first row\n\
of the output is @var{y0}.\n\
\n\
\n\
This is a developmental version of sdirk4().\n\
This function was created from the 1995 fortran source sdirk4.f\n\
which is a Singly Diagonally Implicit Runge Kutta (SDIRK)\n\
ODE integrator specifically designed to solve stiff ODE's.\n\
\n\
Argument descriptions:\n\n\
@var{fcn} is the function that returns the derivatives to be integrated\n\n\
@var{y0} is the vector of initial conditions\n\n\
@var{xspan} is a 2-element vector containing the starting and ending integration times\n\
(i.e. @var{tspan}=[t_start,t_stop]).\n\
@var{rel_tol} & @var{abs_tol} are optional arguments specifying the relative and \n\
absolute error tolerances (default for both is 1e-6)\n\n\
@var{h_initial} is an optional argument specifying the initial stepsize. (default=1e-6)\n\n\
@var{abscissa_trace} is an optional integer argument (0 or 1)that, if 1, will\n\
display every other @var{skip_step}'th integration step and corresponding \n\
independent coordinate value. (default is 'false')\n\n\
@var{skip_step} is an integer and only active when @var{abscissa_trace} is true.\n\n\
note: @var{abscissa_trace} and @var{skip_step} do not affect the variables @var{xout}\n\
or @var{yout}.  They only affect output to the console while sdirk4 is integrating.\n\
\n\n\
example usage:\n\n\
function ydot = statedot(y,x)\n\
\tydot = [-y(1) ; -10000*y(2)];\n\
end\n\n\
\n\n\
x0=[1;2];  tspan=[0 10];  rel_tol=1e-6;  abs_tol=1e-6;\n\
h_initial=1e-3;  abscissa_trace=1;  skip_step=9;\n\
\n\n\
[x,y]=sdirk4('statedot',x0,tspan,rel_tol,abs_tol,h_initial,abscissa_trace,skip_step);\n\n\
axis;  plot(x,y);  pause;  axis([0 .001]);  replot\n\n\
\n\
\n")
{
  octave_value_list retval;

  unwind_protect::begin_frame ("Fsdirk4");

  unwind_protect_int (call_depth);
  call_depth++;

  if (call_depth > 1) {
    ::error( "sdirk4: invalid recursive call") ;
    unwind_protect::run_frame ("Fsdirk4");
    return retval ;
  }

  int nargin = args.length ();

  if (nargin > 2 && nargin < 9 && nargout < 3)
    {
       //fprintf( stderr, "nargin %d nargout %d\n", nargin, nargout ) ;

       // get the user-defined state-dot function
       sdirk4_fcn = extract_function
	   (args(0), "sdirk4:extract_function", "__sdirk4_fcn__",
	    "function xdot = __sdirk4_fcn__ (x, t) xdot = ",
	    "; endfunction");

      if (error_state || ! sdirk4_fcn)
	SDIRK4_ABORT ();

      ColumnVector state (args(1).vector_value ());
      int n = state.capacity ();       // N (int, DIMENSION OF THE SYSTEM)

      if (error_state)
	SDIRK4_ABORT1 ("sdirk4: expecting state vector initial conditions as second argument");

      ColumnVector tspan (args(2).vector_value ());

      if (tspan.capacity() > 2)
        SDIRK4_ABORT1 ("sdirk4: tspan array greater than two elements; tspan should be [t_start,tfinal]");

      t_start = tspan(0);
      t_stop = tspan(1);

      if (error_state)
	SDIRK4_ABORT1 ("sdirk4: expecting tspan vector as second argument, tspan=[t_start,t_final]");

      if ( nargin > 3 ) {
        rel_tol = args(3).double_value();
        if (error_state) {
          unwind_protect::run_frame ("Fsdirk4");
          ::error( "sdirk4: expecting scalar real value (rel_tol) as optional fourth argument");
          return retval ;
        }
      }
      else
        rel_tol = 1e-6;

      if ( nargin > 4 ) {
        abs_tol = args(4).double_value();
        if (error_state) {
          unwind_protect::run_frame ("Fsdirk4");
          ::error( "sdirk4: expecting scalar real value (abs_tol) as optional fourth argument");
          return retval ;
        }
      }
      else
        abs_tol = 1e-6;

      if ( nargin > 5 ) {
        h_initial = args(5).double_value();
        if (error_state) {
          unwind_protect::run_frame ("Fsdirk4");
          ::error( "sdirk4: expecting scalar h_initial as optional fifth argument");
          return retval ;
        }
      }
      else
        h_initial = 0.0;
      //fprintf( stderr, "sdirk4.cc: initial stepsize, h_initial = %f\n",h_initial);

      if ( nargin > 6 ) {
        abscissa_trace = args(6).int_value();
        if (error_state) {
          unwind_protect::run_frame ("Fsdirk4");
          ::error( "sdirk4: expecting scalar abscissa_trace (0 or 1) as optional sixth argument");
          return retval ;
        }
      }
      else
        abscissa_trace = false;

      if ( nargin > 7 ) {
        skip_step = args(7).int_value();
        if (error_state) {
          unwind_protect::run_frame ("Fsdirk4");
          ::error( "sdirk4: expecting scalar integer skip_step (0 to any large integer) as optional seventh argument");
          return retval ;
        }
      }
      else
        skip_step = 0;
      step_counter = 0;

      if (nargin > 8)
	{
	  SDIRK4_ABORT1 ("sdirk4: too many input arguments");
	}

  // setup the call to fortran
  // sdirk4.f specific values
  h = h_initial;        // INITIAL STEP SIZE GUESS, IF H=0.D0, THE CODE PUTS H=1.D-6
  itol = 0;             // ITOL (int, ITOL=0: BOTH RTOL AND ATOL ARE SCALARS)
  ijac = 0;             // IJAC (int, IJAC=0: JACOBIAN IS COMPUTED INTERNALLY)
  mljac = n;            // MLJAC (int, MLJAC=N: JACOBIAN IS A FULL MATRIX)
  mujac = 0;            // MUJAC (int, NEED NOT BE DEFINED IF MLJAC=N.)
  imas = 0;             // IMAS (int, IF IMAS=0, THIS MATRIX IS ASSUMED TO BE THE IDENTITY)
  mlmas = n;            // MLMAS (int, MLMAS=N: THE FULL MATRIX CASE)
  mumas = 0;            // MUMAS (int, NEED NOT BE DEFINED IF MLMAS=N)
  iout = 1;             // IOUT (int, IF IOUT=1, SOLOUT IS CALLED AFTER EVERY SUCCESSFUL STEP)
  lwork = 2*n*n+12*n+7; // LWORK (int, LWORK = 2*N*N+12*N+7)
  liwork = 2*n+4;       // LIWORK (int, 2*N+4)
  lrcont = 5*n+2;       // LRCONT (int, 5*N+2)
  idid = 0;             // IDID (REPORTS ON SUCCESSFULNESS UPON RETURN,
                        //       IDID=1  COMPUTATION SUCCESSFUL,
                        //       IDID=-1 COMPUTATION UNSUCCESSFUL)

  // allocate memory
  double* work = new double[lwork];
  int* iwork = new int[liwork];
  double* statedouble = new double[n];

  // This is the way to convert information from an octave ColumnVector
  // into a matrix of double's.  You have to do this because you have
  // to call fortran functions with double's & not ColumnVectors.
  for (int i = 0; i < n; i++)
    statedouble[i] = state(i);   // <-- this transitions from octave to c++

  for (int i = 0; i < 4; i++)
    iwork[i] = 0; // initialize integrator to default values (remember
                  // fortran indicies are one greater than c++ indicies)

  for (int i = 0; i < 7; i++)
    work[i] = 0; // initialize integrator to default values

  //fprintf( stderr, "sdirk4.cc: debug pt 1.0\n");
  //for (int i = 0; i < n; i++)
  //  fprintf( stderr, "sdirk4.cc: initial conditions: statedouble[%d]=%f\n",i,statedouble[i]);

  // call the fortran function sdirk4_
  F77_XFCN (sdirk4, SDIRK4, (&n, sdirk4_user_function, 
                             &t_start, statedouble, &t_stop, &h,
                             &rel_tol, &abs_tol, &itol,
                             jac,
                             &ijac, &mljac, &mujac,
                             mas,
                             &imas, &mlmas, &mumas,
                             solout,
                             &iout, work, &lwork, iwork,
                             &liwork, &lrcont, &idid));
  // free memory
  delete work;
  delete iwork;
  delete statedouble;

  //fprintf( stderr, "sdirk4.cc: debug pt 2.0\n");

  // this is a test of your ability to call the user-defined
  // state-dot function that returns the vector dy/dx
  //sdirk4_user_function (&n, &t_start, state.fortran_vec(), retval.fortran_vec());

  //for  (int i=0; i<n; i++)
  //  fprintf( stderr, "sdirk4.cc: x[%d]=%f\n",i,state(i));
  //fprintf( stderr, "rel_tol = %f, & abs_tol = %f\n", rel_tol, abs_tol ) ;

  if (! error_state)
	{
          // resize the retval to make space for nargout #1
          retval.resize (1);
          // assign ColumnVector timehistory
          retval(0) = timehistory;
          if ( nargout > 1 )
              // resize the retval to make space for nargout #2
              retval.resize (2);
	      retval(1) = statehistory;
	}
    }
  else
    print_usage ("sdirk4: print_usage, you might try checking the input arguments\n");

  unwind_protect::run_frame ("Fsdirk4");

  return retval;
}





int solout (int* nr, double* xold, double* x,
              double* y, int* n, int* irtn)
{
  // This function gets called at each successful integration step
  // from within sdirk4.  The idea is to collect the time and states
  // at each call to solout() and tack them onto the bottom row
  // of the ColumnVector timehistory and Matrix and statehistory.

  // each time this function gets called, make space for another row
  // to contain the new states at this value of x, or time.
  statehistory.resize(*nr,*n);
  timehistory.resize(*nr);

  //fprintf( stderr, "sdirk4: This is the %d'th call to function 'solout' at t=%f\n",*nr,*x);

  // copy the fortran-style double* into an octave Matrix
  for (int i = 0; i < *n; i++)
    statehistory(*nr-1,i) = y[i];

  // copy the fortran-style double* into an octave ColumnVector
  timehistory(*nr-1) = *x;

  //fprintf( stderr, "sdirk4: solout: timehistory(%d)=%f\n",*nr-1,timehistory(*nr-1));

  if ((abscissa_trace==true) & (step_counter==skip_step))
  {
    fprintf( stderr, "sdirk4: integration step %d at abscissa = %f\n",*nr,timehistory(*nr-1));
    step_counter=0;
  }
  else
    step_counter++;

   return 0;
};

  // dummy function declaration
  int mas (int*, double*, int*)
  {
   // dummy function
   fprintf( stderr, "sdirk4: mas: call to dummy function 'mas'") ;
   return 0;
  };


  // dummy jacobian function declaration
  int jac (int*, double*, double*, double*, int*)
  {
   // dummy function
   fprintf( stderr, "sdirk4: jac: call to dummy function 'jac'") ;
   return 0;
  };