* add support for sparse-matrix Jacobians ==> works!

* change default in param.ccl to use this
* change default in src/include/config.h to default to no longer
  link in LAPACK routines
* update documentation


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

1 files changed, 50 insertions, 52 deletions

diff --git a/src/elliptic/dense_Jacobian.hh b/src/elliptic/dense_Jacobian.hh
index 1e269d5..0c85f47 100644
--- a/src/elliptic/dense_Jacobian.hh
+++ b/src/elliptic/dense_Jacobian.hh
@@ -1,10 +1,12 @@
-// Jacobian.hh -- data structures for the Jacobian matrix
+// dense_Jacobian.hh -- dense-matrix Jacobian 
 // $Header$
 
 //
 #ifdef HAVE_DENSE_JACOBIAN
-// dense_Jacobian_sparsity -- (no-op) accumulate dense Jacobian sparsity info
-// dense_Jacobian_matrix -- Jacobian stored as a dense matrix
+// dense_Jacobian -- Jacobian stored as a dense matrix
+#endif
+#ifdef HAVE_DENSE_JACOBIAN__LAPACK
+// dense_Jacobian__LAPACK -- dense_Jacobian with LAPACK linear solver
 #endif
 //
 
@@ -18,44 +20,30 @@
 
 #ifdef HAVE_DENSE_JACOBIAN
 //
-// This class accumulates the sparsity pattern of a dense-matrix Jacobian.
-//
-// Since a dense-matrix Jacobian doesn't care about the sparsity pattern,
-// the member functions of this class are just no-ops.
+// This class stores the Jacobian as a dense matrix in Fortran (column)
+// order.
 //
-class	dense_Jacobian_sparsity
-	: public Jacobian_sparsity
+class	dense_Jacobian
+	: public Jacobian
 	{
 public:
-	// record the zeroing of the entire matrix
-	void zero_matrix() { }
+	//
+	// routines to access the matrix
+	//
 
-	// record the setting of a matrix element to a specified value
-	void set_element(int II, int JJ, fp value) { }
+	// get a matrix element
+	fp element(int II, int JJ) const
+		{ return matrix_(JJ,II); }
 
-	// record the summing of a value into a matrix element
-	void sum_into_element(int II, int JJ, fp value) { }
+	// dense matrix ==> all elements are explicitly stored
+	bool is_explicitly_stored(int II, int JJ) const
+		{ return true; }
 
-	// constructor, destructor
-	dense_Jacobian_sparsity(patch_system& ps,
-				bool print_msg_flag = false)
-		: Jacobian_sparsity(ps, print_msg_flag)
-		{ }
-	~dense_Jacobian_sparsity() { }
-	};
-#endif	// HAVE_DENSE_JACOBIAN
 
-//******************************************************************************
+	//
+	// routines for setting values in the matrix
+	//
 
-#ifdef HAVE_DENSE_JACOBIAN
-//
-// This class stores the Jacobian as a dense matrix in Fortran (column)
-// order.
-//
-class	dense_Jacobian_matrix
-	: public Jacobian_matrix
-	{
-public:
 	// zero the entire matrix
 	void zero_matrix();
 
@@ -67,16 +55,30 @@ public:
 	void sum_into_element(int II, int JJ, fp value)
 		{ matrix_(JJ,II) += value; }
 
-	// access (get) a matrix element
-	fp element(int II, int JJ)
-		const
-		{ return matrix_(JJ,II); }
+	//
+	// constructor, destructor
+	//
+protected:
+	dense_Jacobian(patch_system& ps,
+		       bool print_msg_flag = false);
+	~dense_Jacobian() { }
+
+protected:
+	// Fortran storage order ==> subscripts are (JJ,II)
+	jtutil::array2d<fp> matrix_;
+	};
+#endif	// HAVE_DENSE_JACOBIAN
 
-	// this is a dense matrix ==> all values are stored explicitly
-	bool is_explicitly_stored(int II, int JJ)
-		const
-		{ return true; }
+//******************************************************************************
 
+#ifdef HAVE_DENSE_JACOBIAN__LAPACK
+//
+// This class defines the linear solver routine using LAPACK.
+//
+class	dense_Jacobian__LAPACK
+	: public dense_Jacobian
+	{
+public:
 	// solve linear system J.x = rhs via LAPACK
 	// ... rhs and x are nominal-grid gridfns
 	// ... overwrites Jacobian matrix with its LU decomposition
@@ -87,20 +89,16 @@ public:
 
 	// constructor, destructor
 public:
-	dense_Jacobian_matrix(patch_system& ps,
-			      dense_Jacobian_sparsity& JS,
-			      bool print_msg_flag = false);
-	~dense_Jacobian_matrix();
+	dense_Jacobian__LAPACK(patch_system& ps,
+			       bool print_msg_flag = false);
+	~dense_Jacobian__LAPACK();
 
 private:
-	// subscripts are (JJ,II) (both 0-origin)
-	jtutil::array2d<fp> matrix_;
-
 	// pivot vector for LAPACK routines
-	integer *pivot_;
+	integer *pivot_;	// size N_rows_
 
 	// work vectors for LAPACK condition number computation
-	integer *iwork_;
-	fp *rwork_;
+	integer *iwork_;	// size N_rows_
+	fp *rwork_;		// size 4*N_rows_
 	};
-#endif	// HAVE_DENSE_JACOBIAN
+#endif	// HAVE_DENSE_JACOBIAN__LAPACK