# Copyright (C) 2001 by Peter J. Gawthrop

###############################################################
## Version control history
###############################################################
## $Header$
## $Log$
## Revision 1.322  2001/08/01 04:06:07  geraint
## Added -i dassl for -cc and -oct.
##
## Revision 1.321  2001/07/27 23:43:34  geraint
## Added -cc to usage options (required for use with xmtt).
##
## Revision 1.320  2001/07/27 23:38:38  geraint
## Changed comment below (# SUMMARY) - fixes xmtt.
##
## Revision 1.319  2001/07/27 23:29:10  geraint

$1_ode2odes_dassl.o	: $1_ode.o	$1_odeo.o	mtt_dassl.o
	@echo "Creating \$@"
	ar -cr \$@ \$^
$1_ode2odes_dassl.oct	: $1_ode.oct	$1_odeo.oct	mtt_dassl.oct
	@echo > /dev/null

mtt_Solver.cc:	mtt_Solver.hh
mtt_AlgebraicSolver.cc: mtt_AlgebraicSolver.hh
$1_ode2odes_${algebraic_solver}.cc:	mtt_Solver.cc mtt_${algebraic_solver}.hh mtt_${algebraic_solver}.cc
$1_ode2odes_${algebraic_solver}.o:	mtt_Solver.o mtt_${algebraic_solver}.o
$1_ode2odes_${algebraic_solver}.cc:	mtt_Solver.cc mtt_AlgebraicSolver.cc mtt_${algebraic_solver}.hh mtt_${algebraic_solver}.cc
$1_ode2odes_${algebraic_solver}.o:	mtt_Solver.o  mtt_AlgebraicSolver.o  mtt_${algebraic_solver}.o
	@echo "Creating $1_ode2odes_${algebraic_solver}.o"
	ar -cr \$@ \$^

#SUMMARY numpar	numerical parameter declaration (m) 
$1_numpar.m:  $1_numpar.txt $1_sympars.txt
	mtt_txt2m $1 numpar

$1_ode2odes.o: $1_ode2odes.cc $1_ode2odes_common.o $1_ode2odes_${integration_method}.o $1_ode2odes_${algebraic_solver}.o
	echo Creating $1_ode2odes.o
	${MTT_CXX} ${MTT_CXXFLAGS} ${MTT_CXXINCS} -c $1_ode2odes.cc -DSTANDALONE

$1_ode2odes.oct: $1_ode2odes.cc $1_ode2odes_common.oct $1_ode2odes_${integration_method}.oct $1_ode2odes_${algebraic_solver}.o
	touch $1_ode2odes.m
	echo Creating $1_ode2odes.oct
	$MKOCTFILE $1_ode2odes.cc mtt_${algebraic_solver}.cc mtt_Solver.cc
	$MKOCTFILE $1_ode2odes.cc mtt_${algebraic_solver}.cc mtt_Solver.cc mtt_AlgebraicSolver.cc

$1_ode2odes.cc: $1_def.r $1_sympars.txt\
		$1_ode2odes_common.m $1_ode2odes_common.cc\
		$1_ode2odes_${integration_method}.m $1_ode2odes_${integration_method}.cc mtt_Solver.cc mtt_${algebraic_solver}.cc mtt_${algebraic_solver}.hh
		$1_ode2odes_${integration_method}.m $1_ode2odes_${integration_method}.cc\
		mtt_Solver.cc mtt_AlgebraicSolver.cc mtt_${algebraic_solver}.cc mtt_${algebraic_solver}.hh
	touch $1_ode2odes.m
	make_ode2odes $1 cc $integration_method $algebraic_solver

#Conversion of m to p to c
#SUMMARY ode	ordinary differential equations (c)
#SUMMARY ode	ordinary differential equations (p)
#SUMMARY state	state declaration (c) 

#! /bin/sh

     ###################################### 
     ##### Model Transformation Tools #####
    ######################################

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.63  2001/08/07 04:39:24  geraint
## Consolidated dassl and residual functions.
##
## Revision 1.62  2001/08/01 22:14:32  geraint
## Bug fix for dassl.
##
## Revision 1.61  2001/08/01 04:06:07  geraint
## Added -i dassl for -cc and -oct.
##
## Revision 1.60  2001/07/16 22:23:00  geraint

inline ColumnVector
mtt_input (ColumnVector &x,
	   ColumnVector &y,
	   const double &t,
	   ColumnVector &par)
{
  static ${algebraic_solver} ae(mtt_ae,MTTNPAR,MTTNU,MTTNX,MTTNY,MTTNYZ);
  static MTT::${algebraic_solver} ae(MTTNPAR,MTTNU,MTTNX,MTTNY,MTTNYZ);
  static ColumnVector u  (MTTNU);

#ifdef STANDALONE
  u = F${sys}_input (x, y, t, par);
#else
  static octave_value_list args, f;
  args (0) = octave_value (x);

#include "mtt_AlgebraicSolver.hh"

ColumnVector
MTT::AlgebraicSolver::solve (const ColumnVector	&x,
			     const ColumnVector	&u,
			     const double	&t,
			     const ColumnVector	&par)
{
  if (_nyz > 0)
    {
      _x = x;
      _uui.insert(u,0);
      _t = t;
      _par = par;
      _ui = ColumnVector(_nyz,1.0);
      Solve();
      _uui.insert(_ui,_nu);
    }
  else
    {
      _uui = u;
    }
  return _uui;
}

ColumnVector
MTT::AlgebraicSolver::eval (const ColumnVector	&ui)
{
  if (_nyz > 0)
    _uui.insert(ui,_nu);
  return mtt_ae(_x,_uui,_t,_par);
}

#ifndef MTT_ALGEBRAICSOLVER
#define MTT_ALGEBRAICSOLVER


#include "mtt_Solver.hh"


namespace MTT
{
  class AlgebraicSolver : public MTT::Solver
  {
  public:

    AlgebraicSolver (const int npar,
		     const int nu,
		     const int nx,
		     const int ny,
		     const int nyz)
      : MTT::Solver (npar,nu,nx,ny,nyz)
    {;}

    ColumnVector
    solve (const ColumnVector	&x,
	   const ColumnVector	&u,
	   const double		&t,
	   const ColumnVector	&par);

    ColumnVector
    eval (const ColumnVector	&ui);

    virtual ~AlgebraicSolver (void) {};

  protected:

    virtual void
    Solve (void) = 0;
  };
}


extern ColumnVector
mtt_ae(ColumnVector	&x,
       ColumnVector	&u,
       const double	&t,
       ColumnVector	&par);


#endif // MTT_ALGEBRAICSOLVER

#ifndef MTT_HJSOLVER
#define MTT_HJSOLVER


#include "mtt_Solver.hh"
#include "mtt_AlgebraicSolver.hh"


namespace MTT
{
class HJ_Solver : public Solver {

  // http://www.netlib.org/opt/hooke.c
  // Hooke and Jeeves solution
  
public:
  class HJ_Solver : public MTT::AlgebraicSolver
  {
    // http://www.netlib.org/opt/hooke.c
    // Hooke and Jeeves solution

  public:
  
  HJ_Solver (sys_ae ae,
	     const int npar,
	     const int nu,
	     const int nx,
	     const int ny,
	     const int nyz)
    : Solver (ae,npar,nu,nx,ny,nyz)
  { static_ptr = this; VARS = nyz; };
    HJ_Solver (const int npar,
	       const int nu,
	       const int nx,
	       const int ny,
	       const int nyz)
      : MTT::AlgebraicSolver (npar,nu,nx,ny,nyz)
    {
      static_ptr = this;
      VARS = nyz;
    }
  
  static double
  f (double tryUi[], int nyz);
    static double
    f (double tryUi[], int nyz);
  
  ~HJ_Solver (void) {};
    ~HJ_Solver (void) {};

protected:
  protected:

  void
  Solve (void);
    void
    Solve (void);

  double
  best_nearby (double    	delta[],
	       double   	point[],
	       double   	prevbest,
	       int      	nvars);
    double
    best_nearby (double    	delta[],
		 double   	point[],
		 double   	prevbest,
		 int      	nvars);
  
  int
  hooke (int	nvars,			  // MTTNYZ
	 double	startpt[],		  // user's initial guess
	 double	endpt[],		  // result
	 double	rho		= 0.05,	  // geometric shrink factor
	 double	epsilon		= 1e-3,	  // end value stepsize
	 int	itermax		= 5000);  // max # iterations 
    int
    hooke (int		nvars,			  // MTTNYZ
	   double	startpt[],		  // user's initial guess
	   double	endpt[],		  // result
	   double	rho		= 0.05,	  // geometric shrink factor
	   double	epsilon		= 1e-3,	  // end value stepsize
	   int		itermax		= 5000);  // max # iterations 

private:
  private:

  int VARS;

public:
    int VARS;
    
  public:

  static HJ_Solver *static_ptr;

};

    static HJ_Solver *static_ptr;
    
  };
}


#endif // MTT_HJSOLVER

#include "mtt_Hybrd_Solver.hh"

// http://www.netlib.org/minpack/hybrd.f
// used by Octave's fsolve
  
Hybrd_Solver *Hybrd_Solver::static_ptr;


MTT::Hybrd_Solver *MTT::Hybrd_Solver::static_ptr;

MTT::Hybrd_Solver (const int npar,
		   const int nu,
		   const int nx,
		   const int ny,
		   const int nyz)
{
  static_ptr = this;
}

ColumnVector
Hybrd_Solver::f_hybrd (const ColumnVector &tryUi)
MTT::Hybrd_Solver::f_hybrd (const ColumnVector &tryUi)
{
  Hybrd_Solver::static_ptr->_yz = Hybrd_Solver::static_ptr->eval(tryUi);
  return Hybrd_Solver::static_ptr->_yz;
  MTT::Hybrd_Solver::static_ptr->_yz = MTT::Hybrd_Solver::static_ptr->eval(tryUi);
  return MTT::Hybrd_Solver::static_ptr->_yz;
}

void
Hybrd_Solver::Solve (void)
MTT::Hybrd_Solver::Solve (void)
{    
  int info;
  static int input_errors;
  static int user_errors;
  static int convergences;
  static int progress_errors;
  static int limit_errors;

		<< "  converge (" << convergences << ") "
		<< "  progress (" << progress_errors << ") "
		<< "  limit (" << limit_errors << ")"
		<< "  unknown (" << unknown_errors << ")"
		<< "  (max error = " << std::abs(eval(_ui).max()) << ")" << std::endl;
    }
}

#ifndef MTT_HYBRDSOLVER
#define MTT_HYBRDSOLVER
#include "mtt_Solver.hh"


#include <octave/NLEqn.h>
#include "mtt_AlgebraicSolver.hh"


namespace MTT
{
class Hybrd_Solver : public Solver {

  // http://www.netlib.org/minpack/hybrd.f
  // used by Octave's fsolve
  
public:
  class Hybrd_Solver : public MTT::AlgebraicSolver
  {
    // http://www.netlib.org/minpack/hybrd.f
    // used by Octave's fsolve
    
  public:

  Hybrd_Solver (sys_ae ae,
    Hybrd_Solver (const int npar,
		const int npar,
		const int nu,
		const int nx,
		const int ny,
		const int nyz)
  : Solver (ae,npar,nu,nx,ny,nyz)
		  const int nu,
		  const int nx,
		  const int ny,
		  const int nyz)
      : MTT::AlgebraicSolver (npar,nu,nx,ny,nyz);
  {
    static_ptr = this;
  }

  static ColumnVector
  f_hybrd (const ColumnVector &tryUi);
    static ColumnVector
    f_hybrd (const ColumnVector &tryUi);

  ~Hybrd_Solver (void) {};
    ~Hybrd_Solver (void) {};

protected:
  protected:

  void
  Solve (void);

public:

  static Hybrd_Solver *static_ptr;

};
    void
    Solve (void);
    
  public:
    
    static Hybrd_Solver *static_ptr;
  };
}


#endif // MTT_HYBRDSOLVER

#include "mtt_Reduce_Solver.hh"


void
Reduce_Solver::Solve (void)
MTT::Reduce_Solver::Solve (void)
{
  //  std::cerr << "Error:"
  //	    << " Symbolic solution of equations failed during model build" << std::endl
  //	    << "       Try using one of the other algebraic solution methods" << std::endl;
  ;
}

ColumnVector
Reduce_Solver::solve (const ColumnVector	&x,
		      const ColumnVector	&u,
		      const double		&t,
		      const ColumnVector	&par)
MTT::Reduce_Solver::solve (const ColumnVector	&x,
			   const ColumnVector	&u,
			   const double		&t,
			   const ColumnVector	&par)
{
  return u;
}

#ifndef MTT_REDUCESOLVER
#define MTT_REDUCESOLVER


#include "mtt_Solver.hh"
#include "mtt_AlgebraicSolver.hh"


namespace MTT
{
class Reduce_Solver : public Solver {

  // Dummy class
  // This will not be used unless the Reduce solver has failed earlier
  // in the model build process
  class Reduce_Solver : public MTT::AlgebraicSolver
  {
    // Dummy class
    // This will not be used unless the Reduce solver has failed earlier
    // in the model build process

public:

  Reduce_Solver (sys_ae ae,
  public:
    
    Reduce_Solver (const int npar,
		 const int npar,
		 const int nu,
		 const int nx,
		 const int ny,
		 const int nyz)
    : Solver (ae,npar,nu,nx,ny,nyz)
  { ; };
		   const int nu,
		   const int nx,
		   const int ny,
		   const int nyz)
      : AlgebraicSolver (npar,nu,nx,ny,nyz)
    {;}
	
  void
  Solve (void);

  ColumnVector
  solve (const ColumnVector	&x,
	 const ColumnVector	&u,
	 const double		&t,
	 const ColumnVector	&par);

  ~Reduce_Solver (void) {};

};
   
    void
    Solve (void);
    
    ColumnVector
    solve (const ColumnVector	&x,
	   const ColumnVector	&u,
	   const double		&t,
	   const ColumnVector	&par);
    
    ~Reduce_Solver (void) {};
  };
}


#endif // MTT_REDUCESOLVER 

#include "mtt_Solver.hh"

Solver::Solver (sys_ae ae,
		const int npar,
		const int nu,
		const int nx,
		const int ny,
		const int nyz)
MTT::Solver::Solver (const int npar,
		     const int nu,
		     const int nx,
		     const int ny,
		     const int nyz)
{
  _ae  = ae;
  _np  = npar; 
  _nu  = nu;
  _nx  = nx;
  _ny  = ny;
  _nyz = nyz;

  _uui = ColumnVector (_nu+_nyz);
};

ColumnVector
Solver::solve (const ColumnVector	&x,
	       const ColumnVector	&u,
	       const double		&t,
	       const ColumnVector	&par)
{
  _x = x;
  _uui.insert(u,0);
  _t = t;
  _par = par;
  _ui = ColumnVector(_nyz,1.0);
  Solve ();
  _uui.insert(_ui,_nu);
  return _uui;
}    

ColumnVector
Solver::eval (const ColumnVector &ui)
{
  _uui.insert(ui,_nu);
  return _ae (_x, _uui, _t, _par);
}

#ifndef MTT_SOLVER
#define MTT_SOLVER


#include <cmath>
#include <cstdlib>
#include <iostream>

#include <octave/oct.h>

class Solver {


namespace MTT
{
protected:
  typedef ColumnVector (*sys_ae) // pointer to F${sys}_ae function
  class Solver
    (ColumnVector &,ColumnVector &,const double &t,ColumnVector &);

public:
  {
  public:

  Solver (sys_ae ae,
	  const int npar,
	  const int nu,
	  const int nx,
	  const int ny,
	  const int nyz);
    Solver (const int npar,
	    const int nu,
	    const int nx,
	    const int ny,
	    const int nyz);

  ColumnVector
  solve (const ColumnVector	&x,
	 const ColumnVector	&u,
	 const double		&t,
	 const ColumnVector	&par);
  
  ColumnVector
  eval (const ColumnVector	&ui);

  virtual ~Solver (void) {};
    virtual ~Solver (void) {};
  
protected:
  
  protected:
    
  virtual void
  Solve (void) = 0;

protected:

  ColumnVector			_x;
  ColumnVector			_uui;
  double			_t;
  ColumnVector			_par;

  ColumnVector  		_ui;
  ColumnVector          	_yz;

  int   			_nu;
  int				_np;
  int				_nx;
  int				_ny;
  int				_nyz;
  
    ColumnVector       		_x;
    ColumnVector	       	_uui;
    double			_t;
    ColumnVector	       	_par;
    
    ColumnVector  		_ui;
    ColumnVector          	_yz;
    
    int   			_nu;
    int				_np;
    int				_nx;
    int				_ny;
    int				_nyz;
  };
}
  sys_ae			_ae;

};

  
#endif // MTT_SOLVER