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

14 files changed, 704 insertions, 162 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_;
 	};
diff --git a/src/gr/AHFinderDirect.hh b/src/gr/AHFinderDirect.hh
index 8e5fc8c..42b77c9 100644
--- a/src/gr/AHFinderDirect.hh
+++ b/src/gr/AHFinderDirect.hh
@@ -69,10 +69,26 @@ extern "C"
 void horizon_function(patch_system& ps,
 		      const struct cactus_grid_info& cgi,
 		      const struct geometry_interpolator_info& gii,
-		      bool Jacobian_flag);
+		      bool Jacobian_flag,
+		      jtutil::norm<fp>& H_norms,
+		      bool msg_flag = true);
 
 // horizon_Jacobian.cc
-Jacobian& create_Jacobian(const patch_system& ps,
+Jacobian& create_Jacobian(patch_system& ps,
 			  const char Jacobian_type[]);
-void horizon_Jacobian(patch_system& ps,
-		      Jacobian& Jac);
+void horizon_Jacobian_SD(patch_system& ps,
+			 Jacobian& Jac);
+void horizon_Jacobian_NP(patch_system& ps,
+			 const struct cactus_grid_info& cgi,
+			 const struct geometry_interpolator_info& ii,
+			 Jacobian& Jac,
+			 fp perturbation_amplitude);
+
+// Newton.cc
+void Newton_solve(patch_system& ps,
+		  const struct cactus_grid_info& cgi,
+		  const struct geometry_interpolator_info& gii,
+		  const char Jacobian_type[],
+		  fp perturbation_amplitude,
+		  int max_Newton_iterations,
+		  fp H_norm_for_convergence);
diff --git a/src/gr/cg.hh b/src/gr/cg.hh
index 9008159..cc0ab31 100644
--- a/src/gr/cg.hh
+++ b/src/gr/cg.hh
@@ -102,6 +102,8 @@
 #define partial_H_wrt_partial_dd_h_22	\
 	p.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_dd_h_22, irho,isigma)
 
+#define Delta_h	p.gridfn(nominal_gfns::gfn__Delta_h, irho,isigma)
+
 //******************************************************************************
 
 //
diff --git a/src/gr/driver.cc b/src/gr/driver.cc
index 546a488..ca978ae 100644
--- a/src/gr/driver.cc
+++ b/src/gr/driver.cc
@@ -3,9 +3,10 @@
 //
 // <<<prototypes for functions local to this file>>>
 // AHFinderDirect_driver - top-level driver
+///
+/// setup_Kerr_horizon - set up Kerr horizon in h (Kerr or Kerr-Schild coords)
 /// setup_ellipsoid - setup up a coordiante ellipsoid in h
-/// setup_Kerr_KerrSchild_horizon - set up Kerr/Kerr-Schild horizon in h
-//
+///
 
 #include <stdio.h>
 #include <assert.h>
@@ -46,12 +47,13 @@ using jtutil::error_exit;
 //
 
 namespace {
+void setup_Kerr_horizon(patch_system& ps,
+			fp x_center, fp y_center, fp z_center,
+			fp m, fp a,
+			bool Kerr_Schild_flag);
 void setup_ellipsoid(patch_system& ps,
 		     fp x_center, fp y_center, fp z_center,
 		     fp x_radius, fp y_radius, fp z_radius);
-void setup_Kerr_KerrSchild_horizon(patch_system& ps,
-				   fp x_center, fp y_center, fp z_center,
-				   fp mass, fp spin);
 	  };
 
 //******************************************************************************
@@ -139,7 +141,6 @@ cgi.K_dd_33_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kzz"));
 //
 // create the patch system and initialize the xyz derivative coefficients
 //
-CCTK_VInfo(CCTK_THORNSTRING, "   creating patch system");
 patch_system ps(origin_x, origin_y, origin_z,
 		patch_system::type_of_name(patch_system_type),
 		N_ghost_points, N_overlap_points, delta_drho_dsigma,
@@ -154,7 +155,7 @@ patch_system ps(origin_x, origin_y, origin_z,
 if	(STRING_EQUAL(initial_guess_method, "read from file"))
    then {
 	CCTK_VInfo(CCTK_THORNSTRING,
-		   "   reading initial guess from \"%s\"",
+		   "reading initial guess from \"%s\"",
 		   initial_guess__read_from_file__file_name);
 	ps.read_ghosted_gridfn(ghosted_gfns::gfn__h,
 			       initial_guess__read_from_file__file_name,
@@ -168,13 +169,22 @@ else if (STRING_EQUAL(initial_guess_method, "ellipsoid"))
 			initial_guess__ellipsoid__x_radius,
 			initial_guess__ellipsoid__y_radius,
 			initial_guess__ellipsoid__z_radius);
+else if (STRING_EQUAL(initial_guess_method, "Kerr/Kerr"))
+   then setup_Kerr_horizon(ps,
+			   initial_guess__Kerr_KerrSchild__x_center,
+			   initial_guess__Kerr_KerrSchild__y_center,
+			   initial_guess__Kerr_KerrSchild__z_center,
+			   initial_guess__Kerr_KerrSchild__mass,
+			   initial_guess__Kerr_KerrSchild__spin,
+			   false);	// use Kerr coords
 else if (STRING_EQUAL(initial_guess_method, "Kerr/Kerr-Schild"))
-   then setup_Kerr_KerrSchild_horizon(ps,
-				      initial_guess__Kerr_KerrSchild__x_center,
-				      initial_guess__Kerr_KerrSchild__y_center,
-				      initial_guess__Kerr_KerrSchild__z_center,
-				      initial_guess__Kerr_KerrSchild__mass,
-				      initial_guess__Kerr_KerrSchild__spin);
+   then setup_Kerr_horizon(ps,
+			   initial_guess__Kerr_KerrSchild__x_center,
+			   initial_guess__Kerr_KerrSchild__y_center,
+			   initial_guess__Kerr_KerrSchild__z_center,
+			   initial_guess__Kerr_KerrSchild__mass,
+			   initial_guess__Kerr_KerrSchild__spin,
+			   true);	// use Kerr-Schild coords
 else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 		   "unknown initial_guess_method=\"%s\"!",
 		   initial_guess_method);			/*NOTREACHED*/
@@ -187,18 +197,49 @@ ps.print_ghosted_gridfn_with_xyz(ghosted_gfns::gfn__h,
 //
 // find the apparent horizon
 //
-if      (STRING_EQUAL(method, "horizon"))
+jtutil::norm<fp> H_norms;
+if      (STRING_EQUAL(method, "horizon function"))
    then {
-	horizon_function(ps, cgi, gii, false);
+	horizon_function(ps, cgi, gii, false, H_norms);
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "   H(h) rms-norm %.2e, infinity-norm %.2e\n",
+		   H_norms.rms_norm(), H_norms.infinity_norm());
 	ps.print_gridfn_with_xyz(nominal_gfns::gfn__H,
 				 true, ghosted_gfns::gfn__h,
 				 "H.dat");
 	}
-else if (STRING_EQUAL(method, "horizon Jacobian"))
+else if (STRING_EQUAL(method, "Jacobian"))
    then {
 	Jacobian& Jac = create_Jacobian(ps, Jacobian_type);
-	horizon_function(ps, cgi, gii, true);
-	horizon_Jacobian(ps, Jac);
+	horizon_function(ps, cgi, gii, true, H_norms);
+	horizon_Jacobian_SD(ps, Jac);
+	print_Jacobian(Jacobian_file_name, Jac);
+	}
+else if (STRING_EQUAL(method, "Jacobian test"))
+   then {
+	Jacobian& SD_Jac = create_Jacobian(ps, Jacobian_type);
+	horizon_function(ps, cgi, gii, true, H_norms);
+	horizon_Jacobian_SD(ps, SD_Jac);
+
+	Jacobian& NP_Jac = create_Jacobian(ps, Jacobian_type);
+	horizon_function(ps, cgi, gii, true, H_norms);
+	horizon_Jacobian_NP(ps, cgi, gii,
+			    NP_Jac,
+			    NP_Jacobian__perturbation_amplitude);
+
+	print_Jacobians(Jacobian_file_name, SD_Jac, NP_Jac);
+	}
+else if (STRING_EQUAL(method, "Newton (NP Jacobian)"))
+   then {
+	Newton_solve(ps, cgi, gii,
+		     Jacobian_type,
+		     NP_Jacobian__perturbation_amplitude,
+		     max_Newton_iterations,
+		     H_norm_for_convergence);
+	ps.print_ghosted_gridfn_with_xyz(ghosted_gfns::gfn__h,
+					 true, ghosted_gfns::gfn__h,
+					 "h.dat",
+					 false);	// no ghost zones
 	}
 else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 		   "unknown method=\"%s\"!",
@@ -206,6 +247,56 @@ else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 }
 
 //******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// This function sets up the horizon of a Kerr black hole in Kerr or
+// Kerr-Schild coordinates, on the nominal grid, in the  h  gridfn.
+//
+// Kerr-Schild coordinates are described in MTW Exercise 33.8, page 903,
+// and the horizon is worked out on page 13.2 of my AHFinderDirect notes.
+//
+// Arguments:
+// [xyz]_center = The position of the Kerr black hole.
+// (m,a) = Describe the Kerr black hole.  Note that my convention has
+//	   a=J/m^2 dimensionless, while MTW take a=J/m=m*(my a).
+// Kerr_Schild_flag = false to use Kerr coordinates,
+//		      true to use Kerr-Schild coordinates
+//
+namespace {
+void setup_Kerr_horizon(patch_system& ps,
+			fp x_center, fp y_center, fp z_center,
+			fp m, fp a,
+			bool Kerr_Schild_flag)
+{
+const char* const name = Kerr_Schild_flag ? "Kerr-Schild" : "Kerr";
+
+CCTK_VInfo(CCTK_THORNSTRING,
+	   "setting up horizon for Kerr in %s coords",
+	   name);
+CCTK_VInfo(CCTK_THORNSTRING,
+	   "   mass=%g, spin=J/m^2=%g, posn=(%g,%g,%g)",
+	   double(m), double(a),
+	   double(x_center), double(y_center), double(z_center));
+
+// horizon in Kerr coordinates is coordinate sphere
+const fp r = m * (1.0 + sqrt(1.0 - a*a));
+
+// horizon in Kerr-Schild coordinates is coordinate ellipsoid
+const fp  z_radius = r;
+const fp xy_radius = Kerr_Schild_flag ? r * sqrt(1.0 + a*a*m*m/(r*r)) : r;
+
+CCTK_VInfo(CCTK_THORNSTRING,
+	   "   horizon is coordinate %s",
+	   Kerr_Schild_flag ? "ellipsoid" : "sphere");
+setup_ellipsoid(ps,
+		x_center, y_center, z_center,
+		xy_radius, xy_radius, z_radius);
+}
+	  }
+
+//******************************************************************************
 
 //
 // This function sets up an ellipsoid in the gridfn h, using the
@@ -234,10 +325,10 @@ void setup_ellipsoid(patch_system& ps,
 		     fp x_radius, fp y_radius, fp z_radius)
 {
 CCTK_VInfo(CCTK_THORNSTRING,
-	   "   h = ellipsoid: center=(%g,%g,%g)",
+	   "setting h = ellipsoid: center=(%g,%g,%g)",
 	   double(x_center), double(y_center), double(z_center));
 CCTK_VInfo(CCTK_THORNSTRING,
-	   "                  radius=(%g,%g,%g)",
+	   "                       radius=(%g,%g,%g)",
 	   double(x_radius), double(y_radius), double(z_radius));
 
 	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
@@ -275,7 +366,7 @@ CCTK_VInfo(CCTK_THORNSTRING,
 		   then r = r_minus;
 		else    CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 				   "\n"
-"   expected exactly one r>0 solution, got 0 or 2!\n"
+"   expected exactly one r>0 solution to quadratic, got 0 or 2!\n"
 "   %s patch (irho,isigma)=(%d,%d) ==> (rho,sigma)=(%g,%g)\n"
 "   direction cosines (xcos,ycos,zcos)=(%g,%g,%g)\n"
 "   ==> r_plus=%g r_minus=%g\n"
@@ -293,62 +384,3 @@ CCTK_VInfo(CCTK_THORNSTRING,
 	}
 }
 	  }
-
-//******************************************************************************
-
-//
-// This function sets up the horizon of a Kerr black hole in Kerr-Schild
-// coordinates, on the nominal grid, in the  h  gridfn.
-//
-// Kerr-Schild coordinates are described in MTW Exercise 33.8, page 903,
-// and the horizon is worked out on page 13 of my AHFinderDirect notes.
-//
-// Arguments:
-// [xyz]_center = The position of the Kerr black hole.
-// (mass,spin) = Describe the Kerr black hole.
-//
-// FIXME FIXME: at present [xyz}_center are ignored.
-//
-namespace {
-void setup_Kerr_KerrSchild_horizon(patch_system& ps,
-				   fp x_center, fp y_center, fp z_center,
-				   fp mass, fp spin)
-{
-CCTK_VInfo(CCTK_THORNSTRING,
-	   "   setting up Kerr/Kerr-Schild horizon");
-CCTK_VInfo(CCTK_THORNSTRING,
-	   "      mass=%g, spin=%g, posn=(%g,%g,%g)",
-	   double(mass), double(spin),
-	   double(x_center), double(y_center), double(z_center));
-
-// horizon in Kerr-Schild is coordinate sphere $r = (1 + \sqrt{1-a^2}) m$
-const fp r = (1.0 + sqrt(1.0 - spin*spin)) * mass;
-
-	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
-	{
-	patch& p = ps.ith_patch(pn);
-
-		for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho)
-		{
-		for (int isigma = p.min_isigma() ;
-		     isigma <= p.max_isigma() ;
-		     ++isigma)
-		{
-		const fp rho   = p.rho_of_irho(irho);
-		const fp sigma = p.sigma_of_isigma(isigma);
-
-		fp theta, phi;
-		p.theta_phi_of_rho_sigma(rho,sigma, theta,phi);
-		fp Kerr_x, Kerr_y, Kerr_z;
-		p.xyz_of_r_rho_sigma(r,rho,sigma, Kerr_x,Kerr_y,Kerr_z);
-
-		const fp KS_x = Kerr_x - mass*spin*sin(theta)*sin(phi);
-		const fp KS_y = Kerr_y + mass*spin*sin(theta)*cos(phi);
-		const fp KS_z = Kerr_z;
-		p.ghosted_gridfn(ghosted_gfns::gfn__h,irho,isigma)
-			= jtutil::hypot3(KS_x, KS_y, KS_z);
-		}
-		}
-	}
-}
-	  }
diff --git a/src/gr/gfn.hh b/src/gr/gfn.hh
index 24c56fe..9f54b5c 100644
--- a/src/gr/gfn.hh
+++ b/src/gr/gfn.hh
@@ -53,6 +53,7 @@ enum	{
 	gfn__partial_H_wrt_partial_dd_h_11,
 	gfn__partial_H_wrt_partial_dd_h_12,
 	gfn__partial_H_wrt_partial_dd_h_22,
-	max_gfn = gfn__partial_H_wrt_partial_dd_h_22	// no comma
+	gfn__Delta_h,
+	max_gfn = gfn__Delta_h		// no comma
 	};
 	  };
diff --git a/src/gr/horizon_Jacobian.cc b/src/gr/horizon_Jacobian.cc
index 7718796..bcaec78 100644
--- a/src/gr/horizon_Jacobian.cc
+++ b/src/gr/horizon_Jacobian.cc
@@ -3,9 +3,11 @@
 //
 // <<<prototypes for functions local to this file>>>
 // create_Jacobian - create a Jacobian matrix of the appropriate type
-// horizon_Jacobian - top-level driver to compute the Jacobian
+//
+// horizon_Jacobian_SD - compute the Jacobian by symbolic differentiation
 /// add_ghost_zone_Jacobian - add ghost zone dependencies to Jacobian
 //
+// horizon_Jacobian_NP - compute the Jacobian by numerical perturbation
 
 #include <stdio.h>
 #include <assert.h>
@@ -53,12 +55,18 @@ void add_ghost_zone_Jacobian(const patch_system& ps,
 	  }
 
 //******************************************************************************
+//******************************************************************************
+//******************************************************************************
 
 //
 // This function decodes the  Jacobian_type  parameter and creates the
 // appropriate type of Jacobian matrix object.
 //
-Jacobian& create_Jacobian(const patch_system& ps,
+// 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& create_Jacobian(patch_system& ps,
 			  const char Jacobian_type[])
 {
 if	(STRING_EQUAL(Jacobian_type, "dense matrix"))
@@ -69,12 +77,14 @@ else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 }
 
 //******************************************************************************
+//******************************************************************************
+//******************************************************************************
 
 //
 // This function computes the Jacobian matrix of the LHS function H(h)
-// from the Jacobian coefficient (angular) gridfns.  The computation is
-// done a Jacobian row at a time, using equation (25) of my 1996 apparent
-// horizon finding paper.
+// by symbolic differentiation from the Jacobian coefficient (angular)
+// gridfns.  The computation is done a Jacobian row at a time, using
+// equation (25) of my 1996 apparent horizon finding paper.
 //
 // Inputs (angular gridfns, on ghosted grid):
 //	h				# shape of trial surface
@@ -84,9 +94,10 @@ else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 // Outputs:
 //	The Jacobian matrix is stored in the Jacobian object Jac.
 //
-void horizon_Jacobian(patch_system& ps, Jacobian& Jac)
+void horizon_Jacobian_SD(patch_system& ps,
+			 Jacobian& Jac)
 {
-CCTK_VInfo(CCTK_THORNSTRING, "   horizon Jacobian");
+CCTK_VInfo(CCTK_THORNSTRING, "   horizon Jacobian_SD");
 
 ps.compute_synchronize_Jacobian();
 Jac.zero();
@@ -189,7 +200,7 @@ Jac.zero();
 	    if (xp.is_in_nominal_grid(xm_irho, xm_isigma))
 	       then Jac(II, xp,xm_irho,xm_isigma) += Jac_rho_sigma;
 	       else add_ghost_zone_Jacobian
-			(ps, xp, 
+			(ps, xp,
 		 xp.corner_ghost_zone_containing_point(m_irho < 0, m_isigma < 0,
 						       xm_irho, xm_isigma),
 			 II, xm_irho, xm_isigma,
@@ -254,3 +265,91 @@ const int ym_ipar_posn = ps.synchronize_Jacobian_y_ipar_posn
 	}
 }
 	  }
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// This function computes the Jacobian matrix of the LHS function H(h)
+// by numerical perturbation of the H(h) function.  The algorithm is as
+// follows:
+//
+// evaluate H(h)
+// array-copy H(h) to scratch array save_H
+//	for each point (y,JJ)
+//	{
+//	const fp save_h_y = h at y;
+//	h at y += perturbation_amplitude;
+//	evaluate H(h)
+//		for each point (x,II)
+//		{
+//		Jac(II,JJ) = (H(II) - save_H(II)) / perturbation_amplitude;
+//		}
+//	h at y = save_h_y;
+//	}
+// array-copy save_H back to H(h)
+//
+void horizon_Jacobian_NP(patch_system& ps,
+			 const struct cactus_grid_info& cgi,
+			 const struct geometry_interpolator_info& ii,
+			 Jacobian& Jac,
+			 fp perturbation_amplitude)
+{
+CCTK_VInfo(CCTK_THORNSTRING, "   horizon Jacobian_NP");
+
+// evaluate H(h)
+jtutil::norm<fp> H_norms;
+horizon_function(ps, cgi, ii, false, H_norms);
+
+// array-copy H(h) to scratch array save_H
+fp *save_H = new fp[Jac.NN()];
+jtutil::array_copy(Jac.NN(), ps.gridfn_data(nominal_gfns::gfn__H), save_H);
+
+	for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn)
+	{
+	patch& yp = ps.ith_patch(ypn);
+	CCTK_VInfo(CCTK_THORNSTRING, "      perturbing in %s patch", yp.name());
+
+	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);
+
+	const fp save_h_y = yp.ghosted_gridfn(ghosted_gfns::gfn__h,
+					      y_irho,y_isigma);
+	yp.ghosted_gridfn(ghosted_gfns::gfn__h, y_irho,y_isigma)
+		+= perturbation_amplitude;
+	H_norms.reset();
+	horizon_function(ps, cgi, ii, false, H_norms, false);
+		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);
+		const fp new_xH = xp.gridfn(nominal_gfns::gfn__H,
+					    x_irho,x_isigma);
+		Jac(II,JJ) = (new_xH - save_H[II]) / perturbation_amplitude;
+		}
+		}
+		}
+	yp.ghosted_gridfn(ghosted_gfns::gfn__h, y_irho,y_isigma)
+		= save_h_y;
+	}
+	}
+   	} 
+
+jtutil::array_copy(Jac.NN(), save_H, ps.gridfn_data(nominal_gfns::gfn__H));
+delete[] save_H;
+}
diff --git a/src/gr/horizon_function.cc b/src/gr/horizon_function.cc
index 82f5f31..9c6cb17 100644
--- a/src/gr/horizon_function.cc
+++ b/src/gr/horizon_function.cc
@@ -46,11 +46,15 @@ using jtutil::error_exit;
 //
 
 namespace {
-void setup_xyz_posns(patch_system& ps);
+void setup_xyz_posns(patch_system& ps, bool msg_flag);
 void interpolate_geometry(patch_system& ps,
 			  const struct cactus_grid_info& cgi,
-			  const struct geometry_interpolator_info& ii);
-void compute_H(patch_system& ps, bool Jacobian_flag);
+			  const struct geometry_interpolator_info& ii,
+			  bool msg_flag);
+void compute_H(patch_system& ps,
+	       bool Jacobian_flag,
+	       jtutil::norm<fp>& H_norms,
+	       bool msg_flag);
 	  }
 
 //******************************************************************************
@@ -92,26 +96,30 @@ void compute_H(patch_system& ps, bool Jacobian_flag);
 // Arguments:
 // Jacobian_flag = true to compute the Jacobian coefficients,
 //		   false to skip this.
+// H_norms = (out) This is set to the gridwise norms of the H(h) function.
 //
 void horizon_function(patch_system& ps,
 		      const struct cactus_grid_info& cgi,
 		      const struct geometry_interpolator_info& ii,
-		      bool Jacobian_flag)
+		      bool Jacobian_flag,
+		      jtutil::norm<fp>& H_norms,
+		      bool msg_flag = true)
 {
-CCTK_VInfo(CCTK_THORNSTRING, "   horizon function");
+if (msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING, "   horizon function");
 
 // fill in values of all ghosted gridfns in ghost zones
 ps.synchronize();
 
 // set up xyz positions of grid points
-setup_xyz_posns(ps);
+setup_xyz_posns(ps, msg_flag);
 
 // interpolate $g_{ij}$, $K_{ij}$ --> $g_{ij}$, $K_{ij}$, $\partial_k g_{ij}$
-interpolate_geometry(ps, cgi, ii);
+interpolate_geometry(ps, cgi, ii, msg_flag);
 
 // compute remaining gridfns --> $H$ and optionally Jacobian coefficients
 // by algebraic ops and angular finite differencing
-compute_H(ps, Jacobian_flag);
+compute_H(ps, Jacobian_flag, H_norms, msg_flag);
 }
 
 //******************************************************************************
@@ -121,9 +129,11 @@ compute_H(ps, Jacobian_flag);
 // in the gridfns global_[xyz].  These will be used by interplate_geometry().
 //
 namespace {
-void setup_xyz_posns(patch_system& ps)
+void setup_xyz_posns(patch_system& ps, bool msg_flag)
 {
-CCTK_VInfo(CCTK_THORNSTRING, "      xyz positions and derivative coefficients");
+if (msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "      xyz positions and derivative coefficients");
 
 	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
 	{
@@ -188,10 +198,12 @@ CCTK_VInfo(CCTK_THORNSTRING, "      xyz positions and derivative coefficients");
 namespace {
 void interpolate_geometry(patch_system& ps,
 			  const struct cactus_grid_info& cgi,
-			  const struct geometry_interpolator_info& gii)
+			  const struct geometry_interpolator_info& gii,
+			  bool msg_flag)
 {
-CCTK_VInfo(CCTK_THORNSTRING,
-	   "      interpolate $g_{ij}$, $K_{ij}$ from Cactus 3-D grid");
+if (msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "      interpolate $g_{ij}$, $K_{ij}$ from Cactus 3-D grid");
 
 int status;
 
@@ -313,8 +325,9 @@ if (first_time)
 	first_time = false;
 
 	// store derivative info in interpolator parameter table
-	CCTK_VInfo(CCTK_THORNSTRING,
-		   "         setting up derivative info for interpolator");
+	if (msg_flag)
+	   then CCTK_VInfo(CCTK_THORNSTRING,
+			   "         setting up derivative info for interpolator");
 
 	status = Util_TableSetIntArray(gii.param_table_handle,
 				       N_OUTPUT_ARRAYS, operand_indices,
@@ -339,9 +352,10 @@ if (first_time)
 			   status);				/*NOTREACHED*/
 	}
 
-CCTK_VInfo(CCTK_THORNSTRING,
-	   "         calling interpolator (%d points)",
-	   N_interp_points);
+if (msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "         calling interpolator (%d points)",
+		   N_interp_points);
 status = CCTK_InterpLocalUniform(N_GRID_DIMS,
 				 gii.operator_handle, gii.param_table_handle,
 				 cgi.coord_origin, cgi.coord_delta,
@@ -421,11 +435,16 @@ if (status < 0)
 // Arguments:
 // Jacobian_flag = true to compute the Jacobian coefficients,
 //		   false to skip this.
+// H_norms = (out) This is set to the gridwise norms of the H(h) function.
 //
 namespace {
-void compute_H(patch_system& ps, bool Jacobian_flag)
+void compute_H(patch_system& ps,
+	       bool Jacobian_flag,
+	       jtutil::norm<fp>& H_norms,
+	       bool msg_flag)
 {
-CCTK_VInfo(CCTK_THORNSTRING, "      computing H(h)");
+if (msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING, "      computing H(h)");
 
 	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
 	{
@@ -504,6 +523,9 @@ CCTK_VInfo(CCTK_THORNSTRING, "      computing H(h)");
 		#include "../gr.cg/horizon_function.c"
 		  }
 
+		// update running norms of H(h) function
+		H_norms.data(H);
+
 		if (Jacobian_flag)
 		   then {
 			// partial_H_wrt_partial_d_h, partial_H_wrt_partial_dd_h
diff --git a/src/gr/make.code.defn b/src/gr/make.code.defn
index 7def0d5..8e9745a 100644
--- a/src/gr/make.code.defn
+++ b/src/gr/make.code.defn
@@ -4,7 +4,8 @@
 # Source files in this directory
 SRCS = driver.cc \
        horizon_function.cc \
-       horizon_Jacobian.cc
+       horizon_Jacobian.cc \
+       Newton.cc
 
 # Subdirectories containing source files
 SUBDIRS =
diff --git a/src/include/config.hh b/src/include/config.hh
index dcbf2c9..6cf3693 100644
--- a/src/include/config.hh
+++ b/src/include/config.hh
@@ -8,7 +8,13 @@
 
 typedef CCTK_REAL fp;
 
+// Fortran 'integer' type
+typedef CCTK_INT integer;
+
 // FIXME: this assumes fp == C 'double'
+//	  (CCTK_REAL_PRECISION_{4,8,16} are helpful, but not quite enough)
+#undef  FP_IS_FLOAT
+#define FP_IS_DOUBLE
 #define FP_SCANF_FORMAT "%lf"
 
 //
diff --git a/src/jtutil/README b/src/jtutil/README
index 102e1b1..1b6ea5e 100644
--- a/src/jtutil/README
+++ b/src/jtutil/README
@@ -13,12 +13,15 @@ linear_map<fp>
 	is a template class (templated on the floating-point type)
 	representing a linear map between a contiguous interval of
 	integers, and floating-point numbers.
+
 cpm_map<fp>
 	is a template class (templated on the floating-point type)
 	representing a mapping from the integers to the integers,
 	of the form  i --> constant +/- i .  (The name abbreviates
 	"Constant Plus or Minus MAP".)
+
 fuzzy<fp>
+	(declarations in "util.hh")
 	is a template class (templated on the floating-point type)
 	providing fuzzy arithmetic, to try to paper over some of the
 	effects of floating-point rounding errors.  For example,
@@ -30,8 +33,20 @@ fuzzy<fp>
 	and have the loop execute as one would naievly expect,
 	even if rounding errors cause the final value of \code{x}
 	to be 0.9999999999999999 or 1.000000000000001 or suchlike.
+
 round<fp>
+	(declarations in "util.hh")
 	is a template class (templated on the floating-point type)
 	to do machine-independent rounding of floating point values
 
+norm<fp>
+	(declarations in "util.hh")
+	is a template class (templated on the floating-point type)
+	to the compute 2-, rms-, and/or infinity-norms of a set of
+	numbers
+
+miscfp.cc
+	(declarations in "util.hh")
+	contains various misc floating-point routines
+
 There are also a number of test drivers in the files  test_*.cc .
diff --git a/src/jtutil/miscfp.cc b/src/jtutil/miscfp.cc
index 5ccda53..4263d2c 100644
--- a/src/jtutil/miscfp.cc
+++ b/src/jtutil/miscfp.cc
@@ -6,12 +6,20 @@
 // jtutil::arctan_xy - 4-quadrant arc tangent
 // jtutil::modulo_reduce - reduce an angle modulo 2*pi radians (360 degrees)
 //
+// jtutil::array_zero - zero an array
+// jtutil::array_copy - copy one array to another
+// jtutil::array_scale - scale an array by a scalar
+//
+// ***** template instantiations *****
+//
 
 #include <math.h>
 #include "stdc.h"
 #include "util.hh"
 
 //******************************************************************************
+//******************************************************************************
+//******************************************************************************
 
 //
 // This function computes the floating point "signum" function (as in APL),
@@ -113,3 +121,78 @@ if (! (fuzzy<double>::GE(xx, xmin) && fuzzy<double>::LE(xx, xmax)) )
 return xx;
 }
 	  }	// namespace jtutil::
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// This function zeros the size-N array dst[].
+//
+namespace jtutil
+	  {
+template <typename fp>
+  void array_zero(int N, fp dst[])
+{
+	for (int i = 0 ; i < N ; ++i)
+	{
+	dst[i] = 0.0;
+	}
+}
+	  }	// namespace jtutil::
+
+//******************************************************************************
+
+//
+// This function copies the size-N array src[] to the size-N array dst[],
+// analogously to the BLAS routines xCOPY().
+//
+namespace jtutil
+	  {
+template <typename fp>
+  void array_copy(int N, const fp src[], fp dst[])
+{
+	for (int i = 0 ; i < N ; ++i)
+	{
+	dst[i] = src[i];
+	}
+}
+	  }	// namespace jtutil::
+
+//******************************************************************************
+
+//
+// This function multiplies the size-N array dst[] by the scalar alpha,
+// analogously to the BLAS routines xSCAL().
+//
+namespace jtutil
+	  {
+template <typename fp>
+  void array_scale(int N, fp alpha, fp dst[])
+{
+	for (int i = 0 ; i < N ; ++i)
+	{
+	dst[i] *= alpha;
+	}
+}
+	  }	// namespace jtutil::
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// ***** template instantiations *****
+//
+
+namespace jtutil
+	  {
+template void array_zero<float >(int, float []);
+template void array_zero<double>(int, double[]);
+
+template void array_copy<float >(int, const float [], float []);
+template void array_copy<double>(int, const double[], double[]);
+
+template void array_scale<float >(int, float , float []);
+template void array_scale<double>(int, double, double[]);
+	  }	// namespace jtutil::
diff --git a/src/jtutil/util.hh b/src/jtutil/util.hh
index 6873e62..89d6bde 100644
--- a/src/jtutil/util.hh
+++ b/src/jtutil/util.hh
@@ -1,13 +1,5 @@
 // util.hh -- stuff for my C++ utility library
 // $Id$
-//
-// jtutil::{how_many_in_range,is_even,is_odd} - misc integer manipulation
-// jtutil::abs - absolute value template
-// jtutil::pow* - raise floating point value to small-integer power
-// jtutil::norm<fp> - compute 2-, rms-, and infinity-norms
-// jtutil::fuzzy::<fp> - fuzzy floating point comparisons and other operations
-// jtutil::round::<fp> - floating point rounding
-//
 
 //
 // prerequisites:
@@ -27,11 +19,15 @@ inline int is_even(int i) { return !(i & 0x1); }
 inline int is_odd (int i) { return  (i & 0x1); }
 
 //
-// absolute value template
+// min/max/absolute value template
 // FIXME: <cmath> puts this in std::, but older g++ versions
 //	  (which don't fully implement namespaces) get confused
 //
 template <typename fp>
+  inline fp min(fp x, fp y) { return (x < y) ? x : y; }
+template <typename fp>
+  inline fp max(fp x, fp y) { return (x > y) ? x : y; }
+template <typename fp>
   inline fp abs(fp x) { return (x > 0.0) ? x : -x; }
 
 //
@@ -66,6 +62,16 @@ double arctan_xy(double x, double y);
 double modulo_reduce(double x, double xmod, double xmin, double xmax);
 
 //
+// simple BLAS-like array arithmetic routines
+//
+template <typename fp>
+  void array_zero(int N, fp dst[]);			// dst[] = 0
+template <typename fp>
+  void array_copy(int N, const fp src[], fp dst[]);	// dst[] = src[]
+template <typename fp>
+  void array_scale(int N, fp alpha, fp dst[]);		// dst[] *= alpha
+
+//
 // more misc math stuff, valid only if <math.h> has been #included
 //
 #ifdef PI	// from "jt/stdc.h"
diff --git a/src/patch/patch_system.cc b/src/patch/patch_system.cc
index dcc54c4..fa0d28d 100644
--- a/src/patch/patch_system.cc
+++ b/src/patch/patch_system.cc
@@ -207,8 +207,11 @@ void patch_system::create_patches(const struct patch_info patch_info_in[],
 				  fp delta_drho_dsigma)
 {
 CCTK_VInfo(CCTK_THORNSTRING,
-	   "constructing %s patch system",
+	   "   constructing %s patch system",
 	   name_of_type(type()));
+CCTK_VInfo(CCTK_THORNSTRING,
+	   "                at origin (%g,%g,%g)",
+	   double(origin_x()), double(origin_y()), double(origin_z()));
 N_grid_points_ = 0;
 ghosted_N_grid_points_ = 0;
 	for (int pn = 0 ; pn < N_patches() ; ++pn)
@@ -216,7 +219,7 @@ ghosted_N_grid_points_ = 0;
 	const struct patch_info& pi = patch_info_in[pn];
 
 	CCTK_VInfo(CCTK_THORNSTRING,
-		   "   constructing %s patch",
+		   "      constructing %s patch",
 		   pi.name);
 
 	struct patch *p;
@@ -322,13 +325,13 @@ const int ghosted_gfn_stride = ghosted_N_grid_points();
 const int N_storage = gfn_stride * N_gridfns_in;
 const int ghosted_N_storage = ghosted_gfn_stride * ghosted_N_gridfns_in;
 
-CCTK_VInfo(CCTK_THORNSTRING, "setting up gridfn storage");
+CCTK_VInfo(CCTK_THORNSTRING, "   setting up gridfn storage");
 CCTK_VInfo(CCTK_THORNSTRING,
-	   "   gfn=[%d,%d] ghosted_gfn=[%d,%d]",
+	   "      gfn=[%d,%d] ghosted_gfn=[%d,%d]",
 	   min_gfn_in, max_gfn_in,
 	   ghosted_min_gfn_in, ghosted_max_gfn_in);
 CCTK_VInfo(CCTK_THORNSTRING,
-	   "   N_grid_points()=%d ghosted_N_grid_points()=%d",
+	   "      N_grid_points()=%d ghosted_N_grid_points()=%d",
 	   N_grid_points(), ghosted_N_grid_points());
 
 // storage arrays for all gridfns
@@ -359,7 +362,7 @@ ghosted_gridfn_storage_ = new fp[ghosted_N_storage];
 	}
 
 CCTK_VInfo(CCTK_THORNSTRING,
-	   "   checking that storage is partitioned properly");
+	   "      checking that storage is partitioned properly");
 
 // check to make sure storage for distinct gridfns
 // forms a partition of the overall storage array
@@ -498,7 +501,7 @@ void patch_system::interlink_full_sphere_patch_system
 	(int N_overlap_points,
 	 int interp_handle, int interp_par_table_handle)
 {
-CCTK_VInfo(CCTK_THORNSTRING, "interlinking full sphere patch system");
+CCTK_VInfo(CCTK_THORNSTRING, "   interlinking full sphere patch system");
 
 patch& pz = ith_patch(patch_number_of_name("+z"));
 patch& px = ith_patch(patch_number_of_name("+x"));
@@ -508,7 +511,7 @@ patch& my = ith_patch(patch_number_of_name("-y"));
 patch& mz = ith_patch(patch_number_of_name("-z"));
 
 // create the ghost zones
-CCTK_VInfo(CCTK_THORNSTRING, "   creating ghost zones");
+CCTK_VInfo(CCTK_THORNSTRING, "      creating ghost zones");
 create_interpatch_ghost_zones(pz, px, N_overlap_points);
 create_interpatch_ghost_zones(pz, py, N_overlap_points);
 create_interpatch_ghost_zones(pz, mx, N_overlap_points);
@@ -523,7 +526,7 @@ create_interpatch_ghost_zones(mz, mx, N_overlap_points);
 create_interpatch_ghost_zones(mz, my, N_overlap_points);
 
 // finish setting up the interpatch ghost zones
-CCTK_VInfo(CCTK_THORNSTRING, "   finishing interpatch setup");
+CCTK_VInfo(CCTK_THORNSTRING, "      finishing interpatch setup");
 finish_interpatch_setup(pz, px,
 			N_overlap_points,
 			interp_handle, interp_par_table_handle);
@@ -574,7 +577,7 @@ void patch_system::interlink_plus_z_hemisphere_patch_system
 	(int N_overlap_points,
 	 int interp_handle, int interp_par_table_handle)
 {
-CCTK_VInfo(CCTK_THORNSTRING, "interlinking +z hemisphere patch system");
+CCTK_VInfo(CCTK_THORNSTRING, "   interlinking +z hemisphere patch system");
 
 patch& pz = ith_patch(patch_number_of_name("+z"));
 patch& px = ith_patch(patch_number_of_name("+x"));
@@ -583,7 +586,7 @@ patch& mx = ith_patch(patch_number_of_name("-x"));
 patch& my = ith_patch(patch_number_of_name("-y"));
 
 // create the ghost zones
-CCTK_VInfo(CCTK_THORNSTRING, "   creating ghost zones");
+CCTK_VInfo(CCTK_THORNSTRING, "      creating ghost zones");
 create_interpatch_ghost_zones(pz, px, N_overlap_points);
 create_interpatch_ghost_zones(pz, py, N_overlap_points);
 create_interpatch_ghost_zones(pz, mx, N_overlap_points);
@@ -598,7 +601,7 @@ mx.create_mirror_symmetry_ghost_zone(mx.min_rho_patch_edge());
 my.create_mirror_symmetry_ghost_zone(my.min_rho_patch_edge());
 
 // finish setting up the interpatch ghost zones
-CCTK_VInfo(CCTK_THORNSTRING, "   finishing interpatch setup");
+CCTK_VInfo(CCTK_THORNSTRING, "      finishing interpatch setup");
 finish_interpatch_setup(pz, px,
 			N_overlap_points,
 			interp_handle, interp_par_table_handle);
@@ -637,7 +640,7 @@ void patch_system::interlink_plus_xy_quadrant_patch_system
 	(int N_overlap_points,
 	 int interp_handle, int interp_par_table_handle)
 {
-CCTK_VInfo(CCTK_THORNSTRING, "interlinking +xy quadrant patch system");
+CCTK_VInfo(CCTK_THORNSTRING, "   interlinking +xy quadrant patch system");
 
 patch& pz = ith_patch(patch_number_of_name("+z"));
 patch& px = ith_patch(patch_number_of_name("+x"));
@@ -645,7 +648,7 @@ patch& py = ith_patch(patch_number_of_name("+y"));
 patch& mz = ith_patch(patch_number_of_name("-z"));
 
 // create the ghost zones
-CCTK_VInfo(CCTK_THORNSTRING, "   creating ghost zones");
+CCTK_VInfo(CCTK_THORNSTRING, "      creating ghost zones");
 create_interpatch_ghost_zones(pz, px, N_overlap_points);
 create_interpatch_ghost_zones(pz, py, N_overlap_points);
 create_interpatch_ghost_zones(px, py, N_overlap_points);
@@ -662,7 +665,7 @@ create_periodic_symmetry_ghost_zones(mz.max_rho_patch_edge(),
 				     true);
 
 // finish setting up the interpatch ghost zones
-CCTK_VInfo(CCTK_THORNSTRING, "   finishing interpatch setup");
+CCTK_VInfo(CCTK_THORNSTRING, "      finishing interpatch setup");
 finish_interpatch_setup(pz, px,
 			N_overlap_points,
 			interp_handle, interp_par_table_handle);
@@ -692,7 +695,7 @@ void patch_system::interlink_plus_xz_quadrant_patch_system
 	(int N_overlap_points,
 	 int interp_handle, int interp_par_table_handle)
 {
-CCTK_VInfo(CCTK_THORNSTRING, "interlinking +xz quadrant patch system");
+CCTK_VInfo(CCTK_THORNSTRING, "   interlinking +xz quadrant patch system");
 
 patch& pz = ith_patch(patch_number_of_name("+z"));
 patch& px = ith_patch(patch_number_of_name("+x"));
@@ -700,7 +703,7 @@ patch& py = ith_patch(patch_number_of_name("+y"));
 patch& my = ith_patch(patch_number_of_name("-y"));
 
 // create the ghost zones
-CCTK_VInfo(CCTK_THORNSTRING, "   creating ghost zones");
+CCTK_VInfo(CCTK_THORNSTRING, "      creating ghost zones");
 create_interpatch_ghost_zones(pz, px, N_overlap_points);
 create_interpatch_ghost_zones(pz, py, N_overlap_points);
 create_interpatch_ghost_zones(pz, my, N_overlap_points);
@@ -717,7 +720,7 @@ create_periodic_symmetry_ghost_zones(py.max_sigma_patch_edge(),
 				     false);
 
 // finish setting up the interpatch ghost zones
-CCTK_VInfo(CCTK_THORNSTRING, "   finishing interpatch setup");
+CCTK_VInfo(CCTK_THORNSTRING, "      finishing interpatch setup");
 finish_interpatch_setup(pz, px,
 			N_overlap_points,
 			interp_handle, interp_par_table_handle);
@@ -747,14 +750,14 @@ void patch_system::interlink_plus_xyz_octant_patch_system
 	(int N_overlap_points,
 	 int interp_handle, int interp_par_table_handle)
 {
-CCTK_VInfo(CCTK_THORNSTRING, "interlinking +xyz octant patch system");
+CCTK_VInfo(CCTK_THORNSTRING, "   interlinking +xyz octant patch system");
 
 patch& pz = ith_patch(patch_number_of_name("+z"));
 patch& px = ith_patch(patch_number_of_name("+x"));
 patch& py = ith_patch(patch_number_of_name("+y"));
 
 // create the ghost zones
-CCTK_VInfo(CCTK_THORNSTRING, "   creating ghost zones");
+CCTK_VInfo(CCTK_THORNSTRING, "      creating ghost zones");
 create_interpatch_ghost_zones(pz, px, N_overlap_points);
 create_interpatch_ghost_zones(pz, py, N_overlap_points);
 create_interpatch_ghost_zones(px, py, N_overlap_points);
@@ -768,7 +771,7 @@ create_periodic_symmetry_ghost_zones(px.min_sigma_patch_edge(),
 				     true);
 
 // finish setting up the interpatch ghost zones
-CCTK_VInfo(CCTK_THORNSTRING, "   finishing interpatch setup");
+CCTK_VInfo(CCTK_THORNSTRING, "      finishing interpatch setup");
 finish_interpatch_setup(pz, px,
 			N_overlap_points,
 			interp_handle, interp_par_table_handle);