add a whole bunch of changes from working at home

--> AHFinderDirect now finds AH correctly for Kerr/offset!!!


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

2 files changed, 274 insertions, 18 deletions

diff --git a/src/elliptic/Jacobian.cc b/src/elliptic/Jacobian.cc
index 9b3e4d8..4d3e411 100644
--- a/src/elliptic/Jacobian.cc
+++ b/src/elliptic/Jacobian.cc
@@ -4,12 +4,22 @@
 //
 // Jacobian::Jacobian
 //
+// print_Jacobian
+// print_Jacobians
+//
 // dense_Jacobian::dense_Jacobian
+// dense_Jacobian::~dense_Jacobian
 // dense_Jacobian::zero
 // dense_Jacobian::zero_row
+// dense_Jacobian::solve_linear_system
 //
 
 #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>
@@ -36,6 +46,42 @@ using jtutil::error_exit;
 #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
+//**************************************
+
+namespace {
+void print_Jacobians_internal(const char file_name[],
+			      const Jacobian& SD_Jac, const Jacobian& NP_Jac,
+			      bool pair_flag);
+	  };
 
 //******************************************************************************
 //******************************************************************************
@@ -44,8 +90,9 @@ using jtutil::error_exit;
 //
 // This function constructs a  Jacobian  object.
 //
-Jacobian::Jacobian(const patch_system& ps)
-	: ps_(ps)
+Jacobian::Jacobian(patch_system& ps)
+	: ps_(ps),
+	  NN_(ps.N_grid_points())
 { }
 
 //******************************************************************************
@@ -53,22 +100,156 @@ Jacobian::Jacobian(const patch_system& ps)
 //******************************************************************************
 
 //
+// This function prints a Jacobian matrix to a named output file.
+//
+void print_Jacobian(const char file_name[], const Jacobian& Jac)
+{
+print_Jacobians_internal(file_name, Jac, Jac, false);
+}
+
+//******************************************************************************
+
+//
+// This function prints a pair of Jacobian matrices (and their difference)
+// to a named output file.
+//
+void print_Jacobians(const char file_name[],
+		     const Jacobian& SD_Jac, const Jacobian& NP_Jac)
+{
+print_Jacobians_internal(file_name, SD_Jac, NP_Jac, true);
+}
+
+//******************************************************************************
+
+//
+// If pair_flag = false, this prints SD_Jac.
+// If pair_flag = true, this prints both Jacobians, and the error in SD_Jac.
+//
+namespace {
+void print_Jacobians_internal(const char file_name[],
+			      const Jacobian& SD_Jac, const Jacobian& NP_Jac,
+			      bool pair_flag)
+{
+const patch_system& ps = SD_Jac.my_patch_system();
+
+FILE *fileptr = fopen(file_name, "w");
+if (fileptr == NULL)
+   then error_exit(ERROR_EXIT,
+"***** dense_Jacobian::print(): can't open output file \"%s\"!",
+		   file_name);					/*NOTREACHED*/
+
+fprintf(fileptr, "# column 1 = x II\n");
+fprintf(fileptr, "# column 2 = x patch number\n");
+fprintf(fileptr, "# column 3 = x irho\n");
+fprintf(fileptr, "# column 4 = x isigma\n");
+fprintf(fileptr, "# column 5 = y JJ\n");
+fprintf(fileptr, "# column 6 = y patch number\n");
+fprintf(fileptr, "# column 7 = y irho\n");
+fprintf(fileptr, "# column 8 = y isigma\n");
+if (pair_flag)
+   then {
+	fprintf(fileptr, "# column 9 = SD_Jac(II,JJ)\n");
+	fprintf(fileptr, "# column 10 = NP_Jac(II,JJ)\n");
+	fprintf(fileptr, "# column 11 = abs error\n");
+	fprintf(fileptr, "# column 12 = rel error\n");
+	}
+   else fprintf(fileptr, "# column 9 = Jac(II,JJ)\n");
+
+	for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn)
+	{
+	patch& xp = ps.ith_patch(xpn);
+
+	for (int x_irho = xp.min_irho() ; x_irho <= xp.max_irho() ; ++x_irho)
+	{
+	for (int x_isigma = xp.min_isigma() ;
+	     x_isigma <= xp.max_isigma() ;
+	     ++x_isigma)
+	{
+	const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho,x_isigma);
+
+		for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn)
+		{
+		patch& yp = ps.ith_patch(ypn);
+
+		for (int y_irho = yp.min_irho() ;
+		     y_irho <= yp.max_irho() ;
+		     ++y_irho)
+		{
+		for (int y_isigma = yp.min_isigma() ;
+		     y_isigma <= yp.max_isigma() ;
+		     ++y_isigma)
+		{
+		const int JJ = ps.gpn_of_patch_irho_isigma(yp, y_irho,y_isigma);
+
+		if (! SD_Jac.is_explicitly_stored(II,JJ))
+		   then continue;			// skip sparse points
+
+		const fp SD = SD_Jac(II,JJ);
+		const fp NP = NP_Jac(II,JJ);
+		const fp abs_error = SD - NP;
+
+		if (pair_flag ? ((SD == 0.0) && (NP == 0.0))
+			      : (SD == 0.0))
+		   then continue;		// skip zero values (if == )
+
+		const fp abs_SD = jtutil::abs(SD);
+		const fp abs_NP = jtutil::abs(NP);
+		const fp rel_error = abs_error / jtutil::max(abs_SD, abs_NP);
+
+		fprintf(fileptr,
+			"%d %d %d %d\t%d %d %d %d\t",
+			II, xpn, x_irho, x_isigma,
+			JJ, ypn, y_irho, y_isigma);
+		if (pair_flag)
+		   then fprintf(fileptr,
+				"%.10g\t%.10g\t%e\t%e\n",
+				double(SD), double(NP),
+				double(abs_error), double(rel_error));
+		   else fprintf(fileptr,
+				"%.10g\n",
+				double(SD));
+		}
+		}
+		}
+	}
+	}
+	}
+
+fclose(fileptr);
+}
+	  };
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
 // This function constructs a  dense_Jacobian  object.
 //
-dense_Jacobian::dense_Jacobian(const patch_system& ps)
+dense_Jacobian::dense_Jacobian(patch_system& ps)
 	: Jacobian(ps),
-	  matrix_(0, ps.N_grid_points(),
-		  0, ps.N_grid_points())
+	  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()
 {
-// FIXME FIXME
+jtutil::array_zero(matrix_.N_array(), matrix_.data_array());
 }
 
 //******************************************************************************
@@ -78,5 +259,48 @@ void dense_Jacobian::zero()
 //
 void dense_Jacobian::zero_row(int II)
 {
-// FIXME FIXME
+	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()
+//
+jtutil::array_copy(NN(), rhs, x);
+
+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*/
 }
diff --git a/src/elliptic/Jacobian.hh b/src/elliptic/Jacobian.hh
index 2c64171..d948f07 100644
--- a/src/elliptic/Jacobian.hh
+++ b/src/elliptic/Jacobian.hh
@@ -4,6 +4,9 @@
 //
 // Jacobian -- abstract base class to describe a Jacobian matrix
 //
+// print_Jacobian - print a Jacobian
+// print_Jacobians - print two Jacobians
+//
 // dense_Jacobian - Jacobian stored as a dense matrix
 //
 
@@ -47,10 +50,6 @@
 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
@@ -72,22 +71,44 @@ public:
 		return operator()(II,JJ);
 		}
 
+	// is a given element explicitly stored, or implicitly 0 via sparsity
+	virtual bool is_explicitly_stored(int II, int JJ) const = 0;
+
 	// zero the whole matrix or a single row
 	virtual void zero()           = 0;
 	virtual void zero_row(int II) = 0;
 
+	// solve linear system J.x = rhs
+	// ... rhs and x are nominal-grid gridfns
+	// ... may modify Jacobian matrix (eg for LU decomposition)
+	virtual void solve_linear_system(int rhs_gfn, int x_gfn) = 0;
+
 	// basic meta-info
-	const patch_system& my_patch_system() const { return ps_; }
+	patch_system& my_patch_system() const { return ps_; }
+	int NN() const { return NN_; }
 
 	// constructor, destructor
-	Jacobian(const patch_system& ps);
+	// ... for analyze-factor-solve sparse-matrix derived classes,
+	//     construction may include the "analyze" operation
+	Jacobian(patch_system& ps);
 	virtual ~Jacobian() { }
 
-private:
-	const patch_system& ps_;
+protected:
+	patch_system& ps_;
+	int NN_;	// ps_.N_grid_points()
 	};
 
 //******************************************************************************
+
+//
+// These functions print 1 or 2 Jacobians to a named output file;
+// in the latter case they also print the error in the SD Jacobian.
+//
+void print_Jacobian(const char file_name[], const Jacobian& Jac);
+void print_Jacobians(const char file_name[],
+		     const Jacobian& SD_Jac, const Jacobian& NP_Jac);
+
+//******************************************************************************
 //******************************************************************************
 //******************************************************************************
 
@@ -105,16 +126,27 @@ public:
 	fp& operator()(int II, int JJ)				// lvalue
 		{ return matrix_(JJ,II); }
 
+	bool is_explicitly_stored(int II, int JJ) const
+		{ return true; }
+
 	// zero the whole matrix or a single row
 	void zero();
 	void zero_row(int II);
 
+	// solve linear system J.x = rhs
+	// ... rhs and x are nominal-grid gridfns
+	// ... may modify Jacobian matrix (eg for LU decomposition)
+	void solve_linear_system(int rhs_gfn, int x_gfn);
+
 	// constructor, destructor
 public:
-	dense_Jacobian(const patch_system& ps);
-	~dense_Jacobian() { }
+	dense_Jacobian(patch_system& ps);
+	~dense_Jacobian();
 
 private:
-	// subscripts are (JJ,II)
+	// subscripts are (JJ,II) (both 0-origin)
 	jtutil::array2d<fp> matrix_;
+
+	// pivot vector for LAPACK routines
+	integer *pivot_;
 	};