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

//
// <<<literal contents of "lapack.hh">>>
//
// decode_Jacobian_type
// new_Jacobian
//
// Jacobian::Jacobian
//
#ifdef HAVE_DENSE_JACOBIAN
// dense_Jacobian::dense_Jacobian
// dense_Jacobian::~dense_Jacobian
// dense_Jacobian::zero_matrix
// dense_Jacobian::zero_row
// dense_Jacobian::solve_linear_system
#endif
//

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

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

#include "stl_vector.hh"

#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"

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

//
// 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.
//#include "lapack.h"
//

//***** begin "lapack.h" contents ******
/* lapack.h -- C/C++ prototypes for (some) BLAS+LAPACK+wrapper routines */
/* $Header$ */

/*
 * prerequisites:
 *	"cctk.h"
 *	"config.hh"
 */

#ifdef __cplusplus
extern "C" {
#endif

/*
 * ***** BLAS *****
 */
integer CCTK_FCALL
  CCTK_FNAME(isamax)(const integer* N, const float SX[], const integer* incx);
integer CCTK_FCALL
  CCTK_FNAME(idamax)(const integer* N, const double DX[], const integer* incx);

/*
 * ***** LAPACK *****
 */
void CCTK_FCALL
  CCTK_FNAME(sgesv)(const integer* N, const integer* NRHS,
		    float A[], const int* LDA,
		    integer IPIV[],
		    float B[], const integer* LDB, integer* info);
void CCTK_FCALL
  CCTK_FNAME(dgesv)(const integer* N, const integer* NRHS,
		    double A[], const int* LDA,
		    integer IPIV[],
		    double B[], const integer* LDB, integer* info);

/*
 * ***** wrappers (for passing character-string args) *****
 */
/* norm_int = 0 for infinity-norm, 1 for 1-norm */
void CCTK_FCALL
  CCTK_FNAME(sgecon_wrapper)(const integer* norm_int,
			     const integer* N,
			     float A[], const integer* LDA,
			     const float* anorm, float* rcond,
			     float WORK[], integer IWORK[],
			     integer* info);
void CCTK_FCALL
  CCTK_FNAME(dgecon_wrapper)(const integer* norm_int,
			     const integer* N,
			     double A[], const integer* LDA,
			     const double* anorm, double* rcond,
			     double WORK[], integer IWORK[],
			     integer* info);

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

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

//
// This function decodes a character string specifying a type (derived class)
// of Jacobian, into an internal enum.
//
enum Jacobian_type
  decode_Jacobian_type(const char Jacobian_type_string[])
{
if	(STRING_EQUAL(Jacobian_type_string, "dense matrix"))
   then {
  #ifdef HAVE_DENSE_JACOBIAN
	return Jacobian_type__dense_matrix;
  #else
	error_exit(ERROR_EXIT,
"\n"
"   decode_Jacobian_type():\n"
"        Jacobian type \"dense matrix\" is not configured in this binary;\n"
"        see \"src/include/config.hh\" for details on what Jacobian types\n"
"        are configured and how to change this\n");		/*NOTREACHED*/
  #endif
	}
else	error_exit(ERROR_EXIT,
"decode_Jacobian_type(): unknown Jacobian_type_string=\"%s\"!\n",
		   Jacobian_type_string);			/*NOTREACHED*/
}

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

//
// This function is an "object factory" for Jacobians: it constructs
// and returns a new-allocated Jacobian object of a specified type.
//
// FIXME: the patch system shouldn't really have to be non-const, but
//	  the Jacobian constructors all require this to allow the
//	  linear solvers to directly update gridfns
//
Jacobian& new_Jacobian(patch_system& ps, enum Jacobian_type Jac_type)
{
switch	(Jac_type)
	{
#ifdef HAVE_DENSE_JACOBIAN
  case Jacobian_type__dense_matrix:
	return *new dense_Jacobian(ps);
#endif
  default:
	error_exit(ERROR_EXIT,
"new_Jacobian(): unknown Jacobian_type=(int)%d!\n",
		   Jac_type);					/*NOTREACHED*/
	}
}

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

//
// This function constructs a  Jacobian  object.
//
Jacobian::Jacobian(patch_system& ps)
	: ps_(ps),
	  NN_(ps.N_grid_points())
{ }

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

#ifdef HAVE_DENSE_JACOBIAN
//
// This function constructs a  dense_Jacobian  object.
//
dense_Jacobian::dense_Jacobian(patch_system& ps)
	: Jacobian(ps),
	  matrix_(0,NN()-1, 0,NN()-1),
	  pivot_(new integer[NN()]),
	  iwork_(new integer[NN()]),
	  rwork_(new fp[4*NN()]) // no comma
{ }
#endif

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

#ifdef HAVE_DENSE_JACOBIAN
//
// THis function destroys a  dense_Jacobian  object.
//
dense_Jacobian::~dense_Jacobian()
{
delete[] iwork_;
delete[] rwork_;
delete[] pivot_;
}
#endif

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

#ifdef HAVE_DENSE_JACOBIAN
//
// This function zeros a  dense_Jacobian  object.
//
void dense_Jacobian::zero_matrix()
{
	for (int JJ = 0 ; JJ < NN() ; ++JJ)
	{
		for (int II = 0 ; II < NN() ; ++II)
		{
		operator()(II,JJ) = 0.0;
		}
	}
}
#endif

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

#ifdef HAVE_DENSE_JACOBIAN
//
// This function zeros a single row of a  dense_Jacobian  object.
//
void dense_Jacobian::zero_row(int II)
{
	for (int JJ = 0 ; JJ < NN() ; ++JJ)
	{
	operator()(II,JJ) = 0.0;
	}
}
#endif

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

#ifdef HAVE_DENSE_JACOBIAN
//
// This function solves the linear system J.x = rhs, with rhs and x
// being nominal-grid gridfns, using LAPACK LU-decomposition routines.
// It returns the estimated infinity-norm condition number of the linear
// system, or 0.0 if the matrix is numerically singular
//
fp dense_Jacobian::solve_linear_system(int rhs_gfn, int x_gfn)
{
const fp *rhs = my_patch_system().gridfn_data(rhs_gfn);
fp       *x   = my_patch_system().gridfn_data(x_gfn);
fp       *A   = matrix_.data_array();

const integer infinity_norm_flag = 0;
const integer stride = 1;
const integer NRHS = 1;
const integer N = NN();
const integer N2 = NN() * NN();

// compute the infinity-norm of the matrix A
// ... max_posn = 1-origin index of A[] element with largest absolute value
#if   defined(FP_IS_FLOAT)
  const integer max_posn = CCTK_FNAME(isamax)(&N2, A, &stride);
#elif defined(FP_IS_DOUBLE)
  const integer max_posn = CCTK_FNAME(idamax)(&N2, A, &stride);
#else
  #error "don't know fp datatype!"
#endif
const fp A_infnorm = jtutil::abs(A[max_posn-1]);

// LU decompose and solve the linear system
//
// ... [sd]gesv() use an "in out" design, where the same argument
//     is used for both rhs and x ==> we must first copy rhs to x
//
	for (int II = 0 ; II < NN() ; ++II)
	{
	x[II] = rhs[II];
	}
integer info;
#if   defined(FP_IS_FLOAT)
  CCTK_FNAME(sgesv)(&N, &NRHS, A, &N, pivot_, x, &N, &info);
#elif defined(FP_IS_DOUBLE)
  CCTK_FNAME(dgesv)(&N, &NRHS, A, &N, pivot_, x, &N, &info);
#else
  #error "don't know fp datatype!"
#endif

if (info < 0)
   then error_exit(ERROR_EXIT,
"\n"
"***** dense_Jacobian::solve_linear_system(rhs_gfn=%d, x_gfn=%d):\n"
"        error return (bad argument) info=%d from [sd]gesv() LAPACK routine!\n"
,
		   rhs_gfn, x_gfn,
		   int(info));					/*NOTREACHED*/

if (info > 0)
   then return 0.0;					// *** ERROR RETURN ***
							// *** (singular matrix)

// estimate infinity-norm condition number
fp rcond;
#if   defined(FP_IS_FLOAT)
  CCTK_FNAME(sgecon_wrapper)(&infinity_norm_flag,
			     &N, A, &N, &A_infnorm, &rcond,
			     rwork_, iwork_, &info);
#elif defined(FP_IS_DOUBLE)
  CCTK_FNAME(dgecon_wrapper)(&infinity_norm_flag,
			     &N, A, &N, &A_infnorm, &rcond,
			     rwork_, iwork_, &info);
#else
  #error "don't know fp datatype!"
#endif
if (rcond == 0.0)
   then return 0.0;					// *** ERROR RETURN ***
							// *** (singular matrix)
return 1.0/rcond;
}
#endif