start adding support for sparse matrix Jacobians (doesn't work yet)

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

3 files changed, 80 insertions, 15 deletions

diff --git a/src/gr/expansion_Jacobian.cc b/src/gr/expansion_Jacobian.cc
index a843e4e..5faf059 100644
--- a/src/gr/expansion_Jacobian.cc
+++ b/src/gr/expansion_Jacobian.cc
@@ -94,7 +94,7 @@ bool expansion_Jacobian(patch_system& ps,
 			const struct cactus_grid_info& cgi,
 			const struct geometry_info& gi,
 			const struct Jacobian_info& Jacobian_info,
-			bool print_msg_flag)
+			bool print_msg_flag /* = false */)
 {
 switch	(Jacobian_info.Jacobian_method)
 	{
@@ -136,8 +136,8 @@ return true;						// *** NORMAL RETURN ***
 
 //
 // This function computes the Jacobian matrix of the expansion Theta(h)
-// by numerical perturbation of the Theta(h) function.  The algorithm is as
-// follows:
+// by numerical perturbation of the Theta(h) function.  The algorithm is
+// as follows:
 //
 // we assume that Theta = Theta(h) has already been evaluated
 // save_Theta = Theta
@@ -161,6 +161,7 @@ return true;						// *** NORMAL RETURN ***
 //
 // Outputs:
 //	The Jacobian matrix is stored in the Jacobian object Jac.
+//	As implied by the above algorithm, it's traversed by columns.
 //
 // Results:
 // This function returns true for a successful computation, or false
@@ -219,7 +220,7 @@ ps.gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta);
 					       x_irho,x_isigma);
 		const fp new_Theta = xp.gridfn(gfns::gfn__Theta,
 					       x_irho,x_isigma);
-		Jac(II,JJ) = (new_Theta - old_Theta) / epsilon;
+		Jac.set_element(II,JJ, (new_Theta - old_Theta) / epsilon);
 		}
 		}
 		}
@@ -238,9 +239,8 @@ return true;						// *** NORMAL RETURN ***
 //
 // This function computes the partial derivative terms in the Jacobian
 // matrix of the expansion Theta(h), 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.
+// the Jacobian coefficient (angular) gridfns.  The Jacobian is traversed
+// by rows, using equation (25) of my 1996 apparent horizon finding paper.
 //
 // Inputs (angular gridfns, on ghosted grid):
 //	h			# shape of trial surface
@@ -316,7 +316,11 @@ ps.compute_synchronize_Jacobian();
 				   * xp.partial_rho_rho_coeff(m_irho);
 	    const fp Jac_sum = Jac_rho + Jac_rho_rho;
 	    if (xp.is_in_nominal_grid(xm_irho, x_isigma))
-	       then Jac(II, xp,xm_irho,x_isigma) += Jac_sum;
+	       then {
+		    const int xm_JJ
+		       = Jac.II_of_patch_irho_isigma(xp,xm_irho,x_isigma);
+		    Jac.sum_into_element(II, xm_JJ, Jac_sum);
+		    }
 	       else add_ghost_zone_Jacobian
 			(ps, Jac,
 			 Jac_sum,
@@ -338,7 +342,11 @@ ps.compute_synchronize_Jacobian();
 				       * xp.partial_sigma_sigma_coeff(m_isigma);
 	    const fp Jac_sum = Jac_sigma + Jac_sigma_sigma;
 	    if (xp.is_in_nominal_grid(x_irho, xm_isigma))
-	       then Jac(II, xp,x_irho,xm_isigma) += Jac_sum;
+	       then {
+		    const int xm_JJ
+		       = Jac.II_of_patch_irho_isigma(xp, x_irho, xm_isigma);
+		    Jac.sum_into_element(II, xm_JJ, Jac_sum);
+		    }
 	       else add_ghost_zone_Jacobian
 			(ps, Jac,
 			 Jac_sum,
@@ -363,7 +371,11 @@ ps.compute_synchronize_Jacobian();
 	       = Jacobian_coeff_rho_sigma
 		 * xp.partial_rho_sigma_coeff(m_irho, m_isigma);
 	    if (xp.is_in_nominal_grid(xm_irho, xm_isigma))
-	       then Jac(II, xp,xm_irho,xm_isigma) += Jac_rho_sigma;
+	       then {
+		    const int xm_JJ
+		       = Jac.II_of_patch_irho_isigma(xp, xm_irho, xm_isigma);
+		    Jac.sum_into_element(II, xm_JJ, Jac_rho_sigma);
+		    }
 	       else {
 		    const ghost_zone& xmgz
 		       = xp.corner_ghost_zone_containing_point
@@ -428,7 +440,7 @@ const patch_edge& ye = ps.synchronize_Jacobian(xmgz,
 					       y_iperp,
 					       y_posn, min_ym, max_ym,
 					       Jacobian_buffer);
-const patch& yp = ye.my_patch();
+patch& yp = ye.my_patch();
 
 // add the Jacobian contributions from the ym points
 	for (int ym = min_ym ; ym <= max_ym ; ++ym)
@@ -436,8 +448,8 @@ const patch& yp = ye.my_patch();
 	const int y_ipar = y_posn + ym;
 	const int y_irho   = ye.  irho_of_iperp_ipar(y_iperp,y_ipar);
 	const int y_isigma = ye.isigma_of_iperp_ipar(y_iperp,y_ipar);
-	const int JJ = ps.gpn_of_patch_irho_isigma(yp, y_irho, y_isigma);
-	Jac(x_II,JJ) += mol*Jacobian_buffer(ym);
+	const int y_JJ = Jac.II_of_patch_irho_isigma(yp, y_irho, y_isigma);
+	Jac.sum_into_element(x_II, y_JJ, mol*Jacobian_buffer(ym));
 	}
 }
 	  }
@@ -498,7 +510,7 @@ if (! expansion(ps, cgi, gi))
 					      irho,isigma);
 		const fp new_Theta = p.gridfn(gfns::gfn__Theta,
 					      irho,isigma);
-		Jac(II,II) += (new_Theta - old_Theta) / epsilon;
+		Jac.sum_into_element(II,II, (new_Theta - old_Theta) / epsilon);
 		}
 		}
 	}
diff --git a/src/gr/gr.hh b/src/gr/gr.hh
index b5e06f7..43df6ae 100644
--- a/src/gr/gr.hh
+++ b/src/gr/gr.hh
@@ -170,7 +170,7 @@ bool expansion_Jacobian(patch_system& ps,
 			const struct cactus_grid_info& cgi,
 			const struct geometry_info& gi,
 			const struct Jacobian_info& Jacobian_info,
-			bool print_msg_flag);
+			bool print_msg_flag = false);
 
 // Schwarzschild_EF.cc
 void Schwarzschild_EF_geometry(patch_system& ps,
@@ -183,3 +183,9 @@ enum geometry_method
 bool geometry_method_is_interp(enum geometry_method geometry_method);
 enum Jacobian_method
   decode_Jacobian_method(const char Jacobian_method_string[]);
+Jacobian_matrix*
+  create_Jacobian_matrix(patch_system& ps,
+			 const struct cactus_grid_info& cgi,
+			 const struct geometry_info& gi,
+			 const struct Jacobian_info& Jac_info,
+			 bool print_msg_flag = false);
diff --git a/src/gr/misc.cc b/src/gr/misc.cc
index 28cdfbd..3951331 100644
--- a/src/gr/misc.cc
+++ b/src/gr/misc.cc
@@ -4,6 +4,7 @@
 // decode_geometry_method - decode the geometry_method parameter
 // geomery_method_is_interp - does  enum geometry_method  specify interpolation?
 // decode_geometry_method - decode the geometry_method parameter
+// create_Jacobian_matrix - create a new Jacobian_matrix object
 //
 
 #include <stdio.h>
@@ -103,3 +104,49 @@ else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 "decode_Jacobian_method(): unknown Jacobian_method_string=\"%s\"!",
 		   Jacobian_method_string);			/*NOTREACHED*/
 }
+
+//******************************************************************************
+
+//
+// This function creates a new-allocated Jacobian matrix and returns a
+// pointer to it.  In more detail, it constructs a Jacobian_sparsity object,
+// calls expansion() and then expansion_Jacobian() to traverse the Jacobian
+// matrix (recording the Jacobian's sparsity pattern in the process), then
+// uses this sparsity pattern to construct the Jacobian matrix itself
+// (and destroys the Jacobian_sparsity object).
+// 
+Jacobian_matrix*
+  create_Jacobian_matrix(patch_system& ps,
+			 const struct cactus_grid_info& cgi,
+			 const struct geometry_info& gi,
+			 const struct Jacobian_info& Jac_info,
+			 bool print_msg_flag /* = false */)
+{
+if (print_msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING,
+                   "traversing Jacobian matrix (recording sparsity pattern)");
+Jacobian_sparsity* JS_ptr
+	= new_Jacobian_sparsity(Jac_info.Jacobian_storage_method,
+				ps,
+				print_msg_flag);
+
+// it's safest to compute expansion() before doing the traversal
+// to make sure all the Jacobian coefficients are set to sane values
+expansion(ps,
+	  cgi, gi,
+	  true);			// yes, compute Jacobian coeffs
+
+// now do the actual traversal (==> records the sparsity pattern)
+expansion_Jacobian(ps, *JS_ptr,
+		   cgi, gi, Jac_info);
+
+if (print_msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "constructing Jacobian matrix");
+Jacobian_matrix* Jac_ptr
+	= new_Jacobian_matrix(Jac_info.Jacobian_storage_method,
+			      ps, *JS_ptr,
+			      print_msg_flag);
+delete JS_ptr;
+return Jac_ptr;
+}