# 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_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 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_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 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);

#ifndef MTT_HJSOLVER
#define MTT_HJSOLVER


#include "mtt_AlgebraicSolver.hh"


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

  public:
  

    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);
  
    ~HJ_Solver (void) {};

  protected:

    void
    Solve (void);

    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 

  private:

    int VARS;
    
  public:

    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

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
MTT::Hybrd_Solver::f_hybrd (const ColumnVector &tryUi)
{
  MTT::Hybrd_Solver::static_ptr->_yz = MTT::Hybrd_Solver::static_ptr->eval(tryUi);
  return MTT::Hybrd_Solver::static_ptr->_yz;
}

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 <octave/NLEqn.h>
#include "mtt_AlgebraicSolver.hh"


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

    Hybrd_Solver (const int npar,

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




    static ColumnVector
    f_hybrd (const ColumnVector &tryUi);

    ~Hybrd_Solver (void) {};

  protected:

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


#endif // MTT_HYBRDSOLVER

#include "mtt_Reduce_Solver.hh"


void
MTT::Reduce_Solver::Solve (void)
{



  ;
}

ColumnVector
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_AlgebraicSolver.hh"


namespace MTT
{
  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 (const int npar,

		   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) {};
  };
}


#endif // MTT_REDUCESOLVER 

#include "mtt_Solver.hh"


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

  _np  = npar; 
  _nu  = nu;
  _nx  = nx;
  _ny  = ny;
  _nyz = nyz;

  _uui = ColumnVector (_nu+_nyz);
};

#ifndef MTT_SOLVER
#define MTT_SOLVER


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

#include <octave/oct.h>


namespace MTT
{

  class Solver

  {
  public:


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










    virtual ~Solver (void) {};
  
  protected:
    





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



  
#endif // MTT_SOLVER