// Jacobian.hh -- data structures for the Jacobian matrix
// $Header$

//
// Jacobian -- abstract base class to describe a Jacobian matrix
#ifdef HAVE_DENSE_JACOBIAN
// dense_Jacobian -- Jacobian stored as a dense matrix
#endif
//
// decode_Jacobian_type - decode string into internal enum
// new_Jacobian - factory method
//

//
// prerequisites:
//	"stdc.h"
//	"config.hh"
//	"../jtutil/util.hh"
//	"../jtutil/array.hh"
//	"../jtutil/cpm_map.hh"
//	"../jtutil/linear_map.hh"
//	"coords.hh"
//	"grid.hh"
//	"fd_grid.hh"
//	"patch.hh"
//	"patch_edge.hh"
//	"patch_interp.hh"
//	"ghost_zone.hh"
//	"patch_system.hh"
//

//******************************************************************************

//
// A Jacobian object stores the (a) Jacobian matrix for a patch system.
// Jacobian is an abstract base class, from which we derive a separate
// concrete class for each type of sparsity pattern (more precisely for
// each type of storage scheme) of the Jacobian matrix.
//
// A row/column index of the Jacobian (denoted II/JJ) is a 0-origin grid
// point number within the patch system.
//
// Note that the APIs here implicitly assume there is only a *single* gridfn
// in the Jacobian computation.  (If we had N gridfns for this, then the
// Jacobian would really be a block matrix with N*N blocks.)  This isn't
// very general, but matches our present use in this thorn.
//

class	Jacobian
	{
public:
	// access a matrix element via row/column indices
	virtual const fp& operator()(int II, int JJ) const = 0;	// rvalue
	virtual       fp& operator()(int II, int JJ)       = 0;	// lvalue

	// access a matrix element via row index and column (patch,irho,isigma)
	const fp& operator()					// rvalue
		(int II,   const patch& JJ_patch, int JJ_irho, int JJ_isigma)
		const
		{
		const int JJ
		   = ps_.gpn_of_patch_irho_isigma(JJ_patch, JJ_irho,JJ_isigma);
		return operator()(II,JJ);
		}
	fp& operator()						// lvalue
		(int II,   const patch& JJ_patch, int JJ_irho, int JJ_isigma)
		{
		const int JJ
		   = ps_.gpn_of_patch_irho_isigma(JJ_patch, JJ_irho,JJ_isigma);
		return operator()(II,JJ);
		}

	// is a given element explicitly stored, or implicitly 0 via sparsity
	virtual bool is_explicitly_stored(int II, int JJ) const = 0;

	// zero the whole matrix or a single row
	virtual void zero_matrix() = 0;
	virtual void zero_row(int II) = 0;

	// solve linear system J.x = rhs
	// ... rhs and x are nominal-grid gridfns
	// ... may modify Jacobian matrix (eg for LU decomposition)
	// ... returns 0.0 if matrix is numerically singular
	//             condition number (> 0.0) if known
	//	       -1.0 if condition number is unknown
	virtual fp solve_linear_system(int rhs_gfn, int x_gfn) = 0;

	// basic meta-info
	patch_system& my_patch_system() const { return ps_; }
	int NN() const { return NN_; }

	// constructor, destructor
	// ... for analyze-factor-solve sparse-matrix derived classes,
	//     construction may include the "analyze" operation
	Jacobian(patch_system& ps);
	virtual ~Jacobian() { }

protected:
	patch_system& ps_;
	int NN_;	// ps_.N_grid_points()
	};

//******************************************************************************

#ifdef HAVE_DENSE_JACOBIAN
//
// This class stores the Jacobian as a dense matrix in Fortran (column)
// order.
//
class	dense_Jacobian
	: public Jacobian
	{
public:
	// access a matrix element via row/column indices
	const fp& operator()(int II, int JJ) const		// rvalue
		{ return matrix_(JJ,II); }
	fp& operator()(int II, int JJ)				// lvalue
		{ return matrix_(JJ,II); }

	bool is_explicitly_stored(int II, int JJ) const
		{ return true; }

	// zero the whole matrix or a single row
	void zero_matrix();
	void zero_row(int II);

	// solve linear system J.x = rhs via LAPACK
	// ... rhs and x are nominal-grid gridfns
	// ... overwrites Jacobian matrix with its LU decomposition
	// ... returns 0.0 if matrix is numerically singular, or
	//             condition number (> 0.0)
	fp solve_linear_system(int rhs_gfn, int x_gfn);

	// constructor, destructor
public:
	dense_Jacobian(patch_system& ps);
	~dense_Jacobian();

private:
	// subscripts are (JJ,II) (both 0-origin)
	jtutil::array2d<fp> matrix_;

	// pivot vector for LAPACK routines
	integer *pivot_;

	// work vectors for LAPACK condition number computation
	integer *iwork_;
	fp *rwork_;
	};
#endif	// HAVE_DENSE_JACOBIAN

//******************************************************************************

enum	Jacobian_type
	{
#ifdef HAVE_DENSE_JACOBIAN
	Jacobian_type__dense_matrix	// no comma on last entry in enum
#endif
	};

// decode string into internal enum
enum Jacobian_type
  decode_Jacobian_type(const char Jacobian_type_string[]);

// construct and return new-allocated Jacobian object of specified type
Jacobian& new_Jacobian(patch_system& ps, enum Jacobian_type Jac_type);