storage for Jacobian matrix which knows ../util/ sparsity structure

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@631 f88db872-0e4f-0410-b76b-b9085cfa78c5

2 files changed, 201 insertions, 0 deletions

diff --git a/src/elliptic/Jacobian.cc b/src/elliptic/Jacobian.cc
new file mode 100644
index 0000000..1bd2b5b
--- /dev/null
+++ b/src/elliptic/Jacobian.cc
@@ -0,0 +1,81 @@
+// Jacobian.cc -- data structures for the Jacobian matrix
+// $Id$
+
+//
+// Jacobian::Jacobian
+//
+// dense_Jacobian::dense_Jacobian
+// dense_Jacobian::zero
+// dense_Jacobian::zero_row
+//
+
+#include <stdio.h>
+#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 "coords.hh"
+#include "grid.hh"
+#include "fd_grid.hh"
+#include "patch.hh"
+#include "patch_edge.hh"
+#include "patch_interp.hh"
+#include "ghost_zone.hh"
+#include "patch_system.hh"
+#include "Jacobian.hh"
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// This function constructs a  Jacobian  object.
+//
+Jacobian::Jacobian(const patch_system& ps)
+	: ps_(ps)
+{ }
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// This function constructs a  dense_Jacobian  object.
+//
+dense_Jacobian::dense_Jacobian(const patch_system& ps)
+	: Jacobian(ps),
+	  matrix_(0, ps.N_grid_points(),
+		  0, ps.N_grid_points())
+{ }
+
+//******************************************************************************
+
+//
+// This function zeros a  dense_Jacobian  object.
+//
+void dense_Jacobian::zero()
+{
+// FIXME FIXME
+}
+
+//******************************************************************************
+
+//
+// This function zeros a single row of a  dense_Jacobian  object.
+//
+void dense_Jacobian::zero_row(int II)
+{
+// FIXME FIXME
+}
diff --git a/src/elliptic/Jacobian.hh b/src/elliptic/Jacobian.hh
new file mode 100644
index 0000000..2c64171
--- /dev/null
+++ b/src/elliptic/Jacobian.hh
@@ -0,0 +1,120 @@
+// Jacobian.hh -- data structures for the Jacobian matrix
+// $Id$
+
+//
+// Jacobian -- abstract base class to describe a Jacobian matrix
+//
+// dense_Jacobian - Jacobian stored as a dense matrix
+//
+
+//
+// prerequisites:
+//	"stdc.h"
+//	"config.hh"
+//	"../jtutil/util.hh"
+//	"../jtutil/array.hh"
+//	"../jtutil/cpm_map.hh"
+//	"../jtutil/linear_map.hh"
+//	"coords.hh"
+//	"grid.hh"
+//	"fd_grid.hh"
+//	"patch.hh"
+//	"patch_edge.hh"
+//	"patch_interp.hh"
+//	"ghost_zone.hh"
+//	"patch_system.hh"
+//
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// A Jacobian object stores the (a) Jacobian matrix for a patch system.
+// Jacobian is an abstract base class, from which we derive a separate
+// concrete class for each type of sparsity pattern (more precisely for
+// each type of storage scheme) of the Jacobian matrix.
+//
+// A row/column index of the Jacobian (denoted II/JJ) is a 0-origin grid
+// point number within the patch system.
+//
+// Note that the APIs here implicitly assume there is only a *single* gridfn
+// in the Jacobian computation.  (If we had N gridfns for this, then the
+// Jacobian would really be a block matrix with N*N blocks.)  This isn't
+// very general, but matches our present use in this thorn.
+//
+
+class	Jacobian
+	{
+public:
+	// basic meta-info
+	patch_system& my_patch_system();
+
+public:
+	// access a matrix element via row/column indices
+	virtual const fp& operator()(int II, int JJ) const = 0;	// rvalue
+	virtual       fp& operator()(int II, int JJ)       = 0;	// lvalue
+
+	// access a matrix element via row index and column (patch,irho,isigma)
+	const fp& operator()					// rvalue
+		(int II,   const patch& JJ_patch, int JJ_irho, int JJ_isigma)
+		const
+		{
+		const int JJ
+		   = ps_.gpn_of_patch_irho_isigma(JJ_patch, JJ_irho,JJ_isigma);
+		return operator()(II,JJ);
+		}
+	fp& operator()						// lvalue
+		(int II,   const patch& JJ_patch, int JJ_irho, int JJ_isigma)
+		{
+		const int JJ
+		   = ps_.gpn_of_patch_irho_isigma(JJ_patch, JJ_irho,JJ_isigma);
+		return operator()(II,JJ);
+		}
+
+	// zero the whole matrix or a single row
+	virtual void zero()           = 0;
+	virtual void zero_row(int II) = 0;
+
+	// basic meta-info
+	const patch_system& my_patch_system() const { return ps_; }
+
+	// constructor, destructor
+	Jacobian(const patch_system& ps);
+	virtual ~Jacobian() { }
+
+private:
+	const patch_system& ps_;
+	};
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// This class stores the Jacobian as a dense matrix in Fortran (column)
+// order.
+//
+class	dense_Jacobian
+	: public Jacobian
+	{
+public:
+	// access a matrix element via row/column indices
+	const fp& operator()(int II, int JJ) const		// rvalue
+		{ return matrix_(JJ,II); }
+	fp& operator()(int II, int JJ)				// lvalue
+		{ return matrix_(JJ,II); }
+
+	// zero the whole matrix or a single row
+	void zero();
+	void zero_row(int II);
+
+	// constructor, destructor
+public:
+	dense_Jacobian(const patch_system& ps);
+	~dense_Jacobian() { }
+
+private:
+	// subscripts are (JJ,II)
+	jtutil::array2d<fp> matrix_;
+	};