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

//
// Jacobian::Jacobian
//
// dense_Jacobian::dense_Jacobian
// dense_Jacobian::~dense_Jacobian
// dense_Jacobian::zero
// dense_Jacobian::zero_row
// dense_Jacobian::solve_linear_system
//
// new_Jacobian
//

#include <stdio.h>
using std::fopen;
using std::printf;
using std::fprintf;
using std::fclose;
using std::FILE;
#include <assert.h>
#include <math.h>
#include <vector>

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

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

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

#include "Jacobian.hh"
//#include "lapack.h"
//**************************************
/* lapack.h -- C/C++ prototypes for (some) LAPACK routines */
/* $Id$ */

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

#ifdef __cplusplus
extern "C" {
#endif

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);

#ifdef __cplusplus
           };	/* extern "C" */
#endif
//**************************************

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

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

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

//
// 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()])
{ }

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

//
// THis function destroys a  dense_Jacobian  object.
//
dense_Jacobian::~dense_Jacobian()
{
delete[] pivot_;
}

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

//
// This function zeros a  dense_Jacobian  object.
//
void dense_Jacobian::zero()
{
	for (int JJ = 0 ; JJ < NN() ; ++JJ)
	{
		for (int II = 0 ; II < NN() ; ++II)
		{
		matrix_(JJ,II) = 0.0;
		}
	}
}

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

//
// 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)
	{
	matrix_(JJ,II) = 0.0;
	}
}

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

//
// This function solves the linear system J.x = rhs, with rhs and x
// being nominal-grid gridfns.  The computation is done using the LAPACK
// [sd]gesv() subroutines; which modify the Jacobian matrix J in-place
// for the LU decomposition.
//
void 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);

//
// [sd]gesv() use an "in out" design, where the same argument is used for
// both rhs and x ==> first copy rhs to x so we can pass that to [sd]gesv()
//
	for (int II = 0 ; II < NN() ; ++II)
	{
	x[II] = rhs[II];
	}

integer N = NN();
integer NRHS = 1;
integer info;

#ifdef FP_IS_FLOAT
  CCTK_FNAME(sgesv)(&N, &NRHS, matrix_.data_array(), &N, pivot_, x, &N, &info);
#endif
#ifdef FP_IS_DOUBLE
  CCTK_FNAME(dgesv)(&N, &NRHS, matrix_.data_array(), &N, pivot_, x, &N, &info);
#endif

if (info != 0)
   then error_exit(ERROR_EXIT,
"\n"
"***** dense_Jacobian::solve_linear_system(rhs_gfn=%d, x_gfn=%d):\n"
"        error return info=%d from [sd]gesv() LAPACK routine!\n"
,
		   rhs_gfn, x_gfn,
		   int(info));					/*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,
		       const char Jacobian_type[])
{
if	(STRING_EQUAL(Jacobian_type, "dense matrix"))
   then return *new dense_Jacobian(ps);
else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
		   "unknown Jacobian_type=\"%s\"!",
		   Jacobian_type);				/*NOTREACHED*/
}