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

//
#ifdef HAVE_ROW_SPARSE_JACOBIAN
// row_sparse_Jacobian -- Jacobian stored as row-oriented sparse matrix
#endif
#ifdef HAVE_ROW_SPARSE_JACOBIAN__ILUCG
// row_sparse_Jacobian__ILUCG -- row_sparse_Jacobian with ILUCG linear solver
#endif
//

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

#undef DEBUG_ROW_SPARSE_JACOBIAN	// defined this for *very* detailed
					// printing of data structures

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

#ifdef HAVE_ROW_SPARSE_JACOBIAN
//
// This class stores the Jacobian as a row-oriented sparse matrix.
//
// This sparse matrix format is widely used by sparse matrix libraries.
// The format is as follows:  The matrix is represented by 3 arrays
// (where () denote Fortran 1-origin subscripting):
//	IA(N_rows+1) = Fortran (1-origin) indices in A of the start
//		       of each row's nonzeros; IA(N_rows+1) = N_nonzeros + 1
//	JA(N_nonzeros) = Fortran (1-origin) column numbers of the
//			 corresponding entries in A
//	A(N_nonzeros) = nonzero matrix elements, ordered by rows;
//			within a row the matrix elements may be in
//			any order
// IA and JA are arrays of  integer  so as to be storage-compatible with
// Fortran.
//
// For example, the matrix
//	[ 11.0  12.0  13.0              16.0 ]
//	[ 21.0  22.0        24.0  25.0       ]
//	[       32.0        34.0             ]
//	[       42.0        44.0        46.0 ]
//	[             53.0        55.0       ]
//	[ 61.0              64.0        66.0 ]
// could be represented by the arrays (Fortran indexing)
//	IA( 1) = 1
//	JA( 1) = 1	JA( 2) = 2	JA( 3) = 3	JA( 4) = 6
//	 A( 1) = 11.0	 A( 2) = 12.0	 A( 3) = 13.0	 A( 4) = 16.0
//	IA( 2) = 5
//	JA( 5) = 1	JA( 6) = 2	JA( 7) = 4	JA( 8) = 5
//	 A( 5) = 21.0	 A( 6) = 22.0	 A( 7) = 24.0	 A( 8) = 25.0
//	IA( 3) = 9
//	JA( 9) = 2	JA(10) = 4
//	 A( 9) = 32.0	 A(10) = 34.0
//	IA( 4) = 11
//	JA(11) = 2	JA(12) = 4	JA(13) = 6
//	 A(11) = 42.0	 A(12) = 44.0	 A(13) = 46.0
//	IA( 5) = 14
//	JA(14) = 3	JA(15) = 5
//	 A(14) = 53.0	 A(15) = 55.0
//	IA( 6) = 16
//	JA(16) = 1	JA(17) = 4	JA(18) = 6
//	 A(16) = 61.0	 A(17) = 64.0	 A(18) = 66.0
//	IA( 7) = 19
//
class	row_sparse_Jacobian
	: public Jacobian
	{
public:
	//
	// routines to access the matrix
	//

	// get a matrix element
	fp element(int II, int JJ) const;

	// is the matrix element (II,JJ) stored explicitly?
	bool is_explicitly_stored(int II, int JJ) const
		{ return find_element(II,JJ) > 0; }


	//
	// routines for setting values in the matrix
	//

	// zero the entire matrix
	void zero_matrix();

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

	// sum a value into a matrix element
	void sum_into_element(int II, int JJ, fp value);


	//
	// constructor, destructor
	//
protected:
	// the constructor only uses ps to get the size of the matrix
	// ... print_msg_flag controls messages about data structure setup
	//     during construction and element-setting process
	row_sparse_Jacobian(patch_system& ps,
			    bool print_msg_flag = false);
	~row_sparse_Jacobian();


	//
	// internal routines
	//
protected:
	// return position (1-origin) of (II,JJ) in A_ and JA_
	// return 0 if not found
	int find_element(int II, int JJ) const;

	// insert new array element (II,JJ) in matrix,
	// starting new row and/or growing arrays if necessary
	// return position (1-origin) in A_ and JA_
	int insert_element(int II, int JJ, fp value);

#ifdef DEBUG_ROW_SPARSE_JACOBIAN
	// print data structure (IA_ and JA_ arrays)
	void print_data_structure() const;
#endif

	// access IA[] and JA[] and A[] using 1-origin Fortran indices
	int& fIA(int fII) const
		{
		assert(fII >= 1);
		assert(fII <= current_N_rows_+1);
		return IA_[csub(fII)];
		}
	int& fJA(int fposn) const
		{
		assert(fposn >= 1);
		assert(fposn <= N_nonzeros_);
		return JA_[csub(fposn)];
		}
	fp& fA(int fposn) const
		{
		assert(fposn >= 1);
		assert(fposn <= N_nonzeros_);
		return A_[csub(fposn)];
		}

protected:
	// actual size, size allocated for JA_ and A_
	int N_nonzeros_, N_nonzeros_allocated_;

	// at present rows 1 <= fii <= current_N_rows_ are valid
	int current_N_rows_;

	// --> new[]-allocated arrays
	// ***** due to constructor initialization list ordering,
	// ***** these must be declared *AFTER* N_nonzeros_allocated_
	integer* IA_;		// size N_rows_+1
				// valid for Fortran indices
				//       1 <= fII <= current_N_rows_+1
	// ... these are grown (reallocated) as necessary by  insert_element()
	// ... using STL vector would be much cleaner here, but alas
	//     there are too many broken compilers/linkers in the world
	//     world don't fully grok vector :( :(
	integer* JA_;		// size N_nonzeros_allocated_
				// valid for Fortran indices
				//       1 <= fposn < IA(current_N_rows_)
	fp* A_;			// ditto
	};
#endif	// HAVE_ROW_SPARSE_JACOBIAN

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

#ifdef HAVE_ROW_SPARSE_JACOBIAN__ILUCG
//
// This class defines the linear solver routine using ILUCG.
//
class row_sparse_Jacobian__ILUCG
	: public row_sparse_Jacobian
	{
public:
	// solve linear system J.x = rhs via ILUCG
	// ... 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,
			       bool print_msg_flag);

	// constructor, destructor
public:
	// the constructor only uses ps to get the size of the matrix
	row_sparse_Jacobian__ILUCG(patch_system& ps,
				   bool print_msg_flag = false);
	~row_sparse_Jacobian__ILUCG();

private:
	bool print_msg_flag_;

	// work vectors for ILUCG subroutine
	// ... allocated by  solve_linear_system()  on first call
	integer* itemp_;
	fp*      rtemp_;
	};
#endif	// HAVE_ROW_SPARSE_JACOBIAN__ILUCG