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

//
#ifdef HAVE_DENSE_JACOBIAN
// dense_Jacobian_sparsity -- (no-op) accumulate dense Jacobian sparsity info
// dense_Jacobian_matrix -- Jacobian stored as a dense matrix
#endif
//

//
// prerequisites:
//	"../patch/patch_system.hh"
//	"Jacobian.hh"
//

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

#ifdef HAVE_DENSE_JACOBIAN
//
// This class accumulates the sparsity pattern of a dense-matrix Jacobian.
//
// Since a dense-matrix Jacobian doesn't care about the sparsity pattern,
// the member functions of this class are just no-ops.
//
class	dense_Jacobian_sparsity
	: public Jacobian_sparsity
	{
public:
	// record the zeroing of the entire matrix
	void zero_matrix() { }

	// record the setting of a matrix element to a specified value
	void set_element(int II, int JJ, fp value) { }

	// record the summing of a value into a matrix element
	void sum_into_element(int II, int JJ, fp value) { }

	// constructor, destructor
	dense_Jacobian_sparsity(patch_system& ps,
				bool print_msg_flag = false)
		: Jacobian_sparsity(ps, print_msg_flag)
		{ }
	~dense_Jacobian_sparsity() { }
	};
#endif	// HAVE_DENSE_JACOBIAN

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

#ifdef HAVE_DENSE_JACOBIAN
//
// This class stores the Jacobian as a dense matrix in Fortran (column)
// order.
//
class	dense_Jacobian_matrix
	: public Jacobian_matrix
	{
public:
	// zero the entire matrix
	void zero_matrix();

	// set a matrix element to a specified value
	void set_element(int II, int JJ, fp value)
		{ matrix_(JJ,II) = value; }

	// sum a value into a matrix element
	void sum_into_element(int II, int JJ, fp value)
		{ matrix_(JJ,II) += value; }

	// access (get) a matrix element
	fp element(int II, int JJ)
		const
		{ return matrix_(JJ,II); }

	// this is a dense matrix ==> all values are stored explicitly
	bool is_explicitly_stored(int II, int JJ)
		const
		{ return true; }

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

	// constructor, destructor
public:
	dense_Jacobian_matrix(patch_system& ps,
			      dense_Jacobian_sparsity& JS,
			      bool print_msg_flag = false);
	~dense_Jacobian_matrix();

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