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

//
// <<<liberal contents of "ilucg.h">>>
//
#ifdef HAVE_ROW_SPARSE_JACOBIAN
// row_sparse_Jacobian::row_sparse_Jacobian
// row_sparse_Jacobian::~row_sparse_Jacobian
// row_sparse_Jacobian::element
// row_sparse_Jacobian::zero_matrix
// row_sparse_Jacobian::set_element
// row_sparse_Jacobian::sum_into_element
/// row_sparse_Jacobian::find_element
/// row_sparse_Jacobian::insert_element
  #ifdef DEBUG_ROW_SPARSE_JACOBIAN
/// row_sparse_Jacobian::print_data_structure
  #endif
#endif
//
#ifdef HAVE_ROW_SPARSE_JACOBIAN__ILUCG
// row_sparse_Jacobian__ILUCG::row_sparse_Jacobian__ILUCG
// row_sparse_Jacobian__ILUCG::~row_sparse_Jacobian__ILUCG
// row_sparse_Jacobian__ILUCG::solve_linear_system
#endif
//

#include <stdio.h>
#include <assert.h>
#include <math.h>
#include <string.h>

#include "util_Table.h"
#include "cctk.h"
#include "cctk_Arguments.h"

#include "config.h"
#include "stdc.h"
#include "../jtutil/util.hh"
#include "../jtutil/array.hh"
#include "../jtutil/cpm_map.hh"
#include "../jtutil/linear_map.hh"
using jtutil::error_exit;

#include "../patch/coords.hh"
#include "../patch/grid.hh"
#include "../patch/fd_grid.hh"
#include "../patch/patch.hh"
#include "../patch/patch_edge.hh"
#include "../patch/patch_interp.hh"
#include "../patch/ghost_zone.hh"
#include "../patch/patch_system.hh"

#include "Jacobian.hh"
#include "row_sparse_Jacobian.hh"

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

//
// FIXME:
//   Cactus's CCTK_FCALL() isn't expanded in .h files (this is a bug),
//   so we include the contents of "lapack.h" inline here.  This is a
//   kludge! :( :(
//#include "ilucg.h"
//

//***** begin "ilucg.h" contents ******
/* ilucg.h -- C/C++ prototypes for [sd]ilucg() routines */
/* $Header$ */

/*
 * prerequisites:
 *	"cctk.h"
 *	"config.h"	// for "integer" = Fortran integer
 */

#ifdef __cplusplus
extern "C"
       {
#endif

/*
 * ***** ILUCG *****
 */
integer CCTK_FCALL
  CCTK_FNAME(silucg)(const integer* N,
		     const integer ia[], const integer ja[], const float a[],
		     const float b[], float x[],
		     integer itemp[], float rtemp[],
		     const float* eps, const integer* iter,
		     integer* istatus);
integer CCTK_FCALL
  CCTK_FNAME(dilucg)(const integer* N,
		     const integer ia[], const integer ja[], const double a[],
		     const double b[], double x[],
		     integer itemp[], double rtemp[],
		     const double* eps, const integer* iter,
		     integer* istatus);

/*
 * ***** wrappers (for returning logical results) *****
 */
void CCTK_FCALL
  CCTK_FNAME(silucg_wrapper)
		    (const integer* N,
		     const integer ia[], const integer ja[], const float a[],
		     const float b[], float x[],
		     integer itemp[], float rtemp[],
		     const float* eps, const integer* iter,
		     integer* istatus, integer* ierror);
void CCTK_FCALL
  CCTK_FNAME(dilucg_wrapper)
		    (const integer* N,
		     const integer ia[], const integer ja[], const double a[],
		     const double b[], double x[],
		     integer itemp[], double rtemp[],
		     const double* eps, const integer* iter,
		     integer* istatus, integer* ierror);

#ifdef __cplusplus
       }	/* extern "C" */
#endif
//***** end "ilucg.h" contents ******

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

#ifdef HAVE_ROW_SPARSE_JACOBIAN
//
// This function constructs a  row_sparse_Jacobian  object.
//
row_sparse_Jacobian::row_sparse_Jacobian(patch_system& ps,
					 bool print_msg_flag /* = false */)
	: Jacobian(ps),
	  N_nonzeros_(0),
	  N_nonzeros_allocated_(FD_GRID__MOL_AREA * N_rows_),
				// initial guess, insert_element()
				// will grow if/when necessary
	  current_N_rows_(0),
	  IA_(new integer[N_rows_+1]),
	  JA_(new integer[N_nonzeros_allocated_]),
	   A_(new fp     [N_nonzeros_allocated_])
{
if (print_msg_flag)
   then CCTK_VInfo(CCTK_THORNSTRING,
		   "   row sparse matrix Jacobian (%d rows)",
		   N_rows_);

zero_matrix();
}
#endif	// HAVE_ROW_SPARSE_JACOBIAN

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

#ifdef HAVE_ROW_SPARSE_JACOBIAN
//
// This function destroys a  row_sparse_Jacobian  object.
//
row_sparse_Jacobian::~row_sparse_Jacobian()
{
delete[]  A_;
delete[] JA_;
delete[] IA_;
}
#endif	// HAVE_ROW_SPARSE_JACOBIAN

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

#ifdef HAVE_ROW_SPARSE_JACOBIAN
//
// This function gets a matrix element.
//
fp row_sparse_Jacobian::element(int II, int JJ)
	const
{
const int fposn = find_element(II,JJ);
return (fposn > 0) ? fA(fposn) : 0.0;
}
#endif	// HAVE_ROW_SPARSE_JACOBIAN

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

#ifdef HAVE_ROW_SPARSE_JACOBIAN
//
// This function zeros a  row_sparse_Jacobian  object.
//
void row_sparse_Jacobian::zero_matrix()
{
#ifdef DEBUG_ROW_SPARSE_JACOBIAN
printf("row_sparse_Jacobian::zero_matrix()\n");
#endif
N_nonzeros_ = 0;
current_N_rows_= 0;
fIA(1) = 1;
}
#endif	// HAVE_ROW_SPARSE_JACOBIAN

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

#ifdef HAVE_ROW_SPARSE_JACOBIAN
//
// This function sets a matrix element to a specified value.
//
void row_sparse_Jacobian::set_element(int II, int JJ, fp value)
{
const int fposn = find_element(II,JJ);
if (fposn > 0)
   then fA(fposn) = value;
   else insert_element(II,JJ, value);
}
#endif	// HAVE_ROW_SPARSE_JACOBIAN

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

#ifdef HAVE_ROW_SPARSE_JACOBIAN
//
// This function sums a value into a matrix element.
//
void row_sparse_Jacobian::sum_into_element(int II, int JJ, fp value)
{
const int fposn = find_element(II,JJ);
if (fposn > 0)
   then fA(fposn) += value;
   else insert_element(II,JJ, value);
}
#endif	// HAVE_ROW_SPARSE_JACOBIAN

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

#ifdef HAVE_ROW_SPARSE_JACOBIAN
//
// This function searches our data structures for element (II,JJ).
// If found, it returns the position (1-origin) in A_ and JA_.
// If not found, it returns 0.
//
int row_sparse_Jacobian::find_element(int II, int JJ)
	const
{
const int fII = fsub(II);
const int fJJ = fsub(JJ);

if (fII > current_N_rows_)
   then return 0;				// this row not defined yet

const int fstart = fIA(fII);
const int fstop  = fIA(fII + 1);
	for (int fposn = fstart ; fposn < fstop ; ++fposn)
	{
	if (fJA(fposn) == fJJ)
	   then return fposn;				// found
	}

return 0;						// not found
}
#endif	// HAVE_ROW_SPARSE_JACOBIAN

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

#ifdef HAVE_ROW_SPARSE_JACOBIAN
//
// This function inserts a new element in the matrix, starting a new
// row and/or growing the arrays if necessary.  It should only be called
// if JJ isn't already in this row.
//
int row_sparse_Jacobian::insert_element(int II, int JJ, fp value)
{
const int fII = fsub(II);
const int fJJ = fsub(JJ);

#ifdef DEBUG_ROW_SPARSE_JACOBIAN
printf(
  "row_sparse_Jacobian::insert_element(): fII=%d fJJ=%d\n",
       fII, fJJ);
#endif

if (fII != current_N_rows_)
   then {
      #ifdef DEBUG_ROW_SPARSE_JACOBIAN
	printf("   starting new row\n");
      #endif
	assert(fII == current_N_rows_+1);
	assert(current_N_rows_ < N_rows_);
	++current_N_rows_;
	fIA(current_N_rows_+1) = fIA(current_N_rows_);
      #ifdef DEBUG_ROW_SPARSE_JACOBIAN
	print_data_structure();
      #endif
	}

if (fIA(fII+1) > N_nonzeros_allocated_)
   then {
      #ifdef DEBUG_ROW_SPARSE_JACOBIAN
	printf("   growing arrays from N_nonzeros_allocated_=%d",
	       N_nonzeros_allocated_);
      #endif
	N_nonzeros_allocated_ += (N_nonzeros_allocated_ >> 1);
      #ifdef DEBUG_ROW_SPARSE_JACOBIAN
	printf(" to %d\n", N_nonzeros_allocated_);
      #endif
	integer* const new_JA = new integer[N_nonzeros_allocated_];
	fp     * const  new_A = new fp     [N_nonzeros_allocated_];
		for (int posn = 0 ; posn < N_nonzeros_ ; ++posn)
		{
		new_JA[posn] = JA_[posn];
		 new_A[posn] =  A_[posn];
		}
	delete[]  A_;
	delete[] JA_;
	JA_ = new_JA;
	 A_ =  new_A;
	}

++N_nonzeros_;
const int fposn = fIA(fII+1)++;
fJA(fposn) = fJJ;
 fA(fposn) = value;
#ifdef DEBUG_ROW_SPARSE_JACOBIAN
printf("   storing at fposn=%d (new N_nonzeros_=%d)\n", fposn, N_nonzeros_);
#endif
return fposn;
}
#endif	// HAVE_ROW_SPARSE_JACOBIAN

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

#ifdef HAVE_ROW_SPARSE_JACOBIAN
#ifdef DEBUG_ROW_SPARSE_JACOBIAN
//
// This function prints the sparse-matrix data structure.
//
void row_sparse_Jacobian::print_data_structure()
	const
{
printf("--- begin Jacobian printout\n");
	for (int fII = 1 ; fII <= current_N_rows_ ; ++fII)
	{
	printf("fII=%d", fII);
	const int fstart = fIA(fII);
	const int fstop  = fIA(fII + 1);
	printf("\tfposn=[%d,%d):", fstart, fstop);
	printf("\tfJJ =");
		for (int fposn = fstart ; fposn < fstop ; ++fposn)
		{
		printf(" %d", fJA(fposn));
		}
	printf("\n");
	}
printf("--- end Jacobian printout\n");
}
#endif	// DEBUG_ROW_SPARSE_JACOBIAN
#endif	// HAVE_ROW_SPARSE_JACOBIAN

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

#ifdef HAVE_ROW_SPARSE_JACOBIAN__ILUCG
//
// This function constructs a  row_sparse_Jacobian__ILUCG  object.
//
row_sparse_Jacobian__ILUCG::row_sparse_Jacobian__ILUCG
	(patch_system& ps,
	 bool print_msg_flag /* = false */)
	: row_sparse_Jacobian(ps, print_msg_flag),
	  print_msg_flag_(print_msg_flag),
	  itemp_(NULL),
	  rtemp_(NULL)
{
if (print_msg_flag)
   then CCTK_VInfo(CCTK_THORNSTRING,
		   "      ILUCG linear-equations solver");
}
#endif	// HAVE_ROW_SPARSE_JACOBIAN__ILUCG

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

#ifdef HAVE_ROW_SPARSE_JACOBIAN__ILUCG
//
// This function destroys a  row_sparse_Jacobian__ILUCG  object.
//
row_sparse_Jacobian__ILUCG::~row_sparse_Jacobian__ILUCG()
{
delete[] rtemp_;
delete[] itemp_;
}
#endif	// HAVE_ROW_SPARSE_JACOBIAN__ILUCG

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

#ifdef HAVE_ROW_SPARSE_JACOBIAN__ILUCG
//
// This function solves the linear system J.x = rhs, with rhs and x
// being nominal-grid gridfns, using an ILUCG subroutine originally
// provided to me in 1985 (!) by Tom Nicol of the UBC Computing Center.
// It returns -1.0 (no condition number estimate).
//
fp row_sparse_Jacobian__ILUCG::solve_linear_system(int rhs_gfn, int x_gfn,
						   bool print_msg_flag)
{
assert(current_N_rows_ == N_rows_);

//
// if this is our first call, allocate the scratch arrays
//
if (itemp_ == NULL)
   then {
	if (print_msg_flag)
	   then {
		CCTK_VInfo(CCTK_THORNSTRING,
			   "row_sparse_Jacobian__ILUCG::solve_linear_system()");
		CCTK_VInfo(CCTK_THORNSTRING,
			   "   N_rows_=%d N_nonzeros_=%d",
			   N_rows_, N_nonzeros_);
		CCTK_VInfo(CCTK_THORNSTRING,
			   "   N_nonzeros_allocated_=%d",
			   N_nonzeros_allocated_);
		}
	itemp_ = new integer[3*N_rows_ + 3*N_nonzeros_ + 2];
	rtemp_ = new fp     [4*N_rows_ +   N_nonzeros_    ];
	}

//
// set up the ILUCG subroutine
//

// initial guess = all zeros
fp *x = ps_.gridfn_data(x_gfn);
	for (int II = 0 ; II < N_rows_ ; ++II)
	{
	x[II] = 0.0;
	}

const integer N = N_rows_;
const fp *rhs = ps_.gridfn_data(rhs_gfn);
const fp eps = -1000.0;			// solve to 1000 * machine epsilon
const integer max_iterations = 0;	// no limit on iterations
integer istatus, ierror;

// the actual linear solution
#if   defined(FP_IS_FLOAT)
  CCTK_FNAME(silucg_wrapper)(&N,
			     IA_, JA_, A_,
			     rhs, x,
			     itemp_, rtemp_,
			     &eps, &max_iterations,
			     &istatus, &ierror);
#elif defined(FP_IS_DOUBLE)
  CCTK_FNAME(dilucg_wrapper)(&N,
			     IA_, JA_, A_,
			     rhs, x,
			     itemp_, rtemp_,
			     &eps, &max_iterations,
			     &istatus, &ierror);
#else
  #error "don't know fp datatype!"
#endif

if (ierror != 0)
   then CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
"***** row_sparse_Jacobian__ILUCG::solve_linear_system(rhs_gfn=%d, x_gfn=%d):\n"
"        error return istatus=%d from [sd]ilucg() routine!"
		   ,
		   rhs_gfn, x_gfn,
		   int(istatus));				/*NOTREACHED*/
if (print_msg_flag)
   then CCTK_VInfo(CCTK_THORNSTRING,
		   "   solution found in %d CG iterations",
		   istatus);

return -1.0;			// no condition number estimate available
}
#endif	// HAVE_ROW_SPARSE_JACOBIAN__ILUCG