* 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, 351 insertions, 104 deletions

diff --git a/src/elliptic/row_sparse_Jacobian.cc b/src/elliptic/row_sparse_Jacobian.cc
index 131655c..7d0be44 100644
--- a/src/elliptic/row_sparse_Jacobian.cc
+++ b/src/elliptic/row_sparse_Jacobian.cc
@@ -2,22 +2,28 @@
 // $Header$
 
 //
+// <<<liberal contents of "ilucg.h">>>
+//
 #ifdef HAVE_ROW_SPARSE_JACOBIAN
-// row_sparse_Jacobian_sparsity::row_sparse_Jacobian_sparsity
-// row_sparse_Jacobian_sparsity::~row_sparse_Jacobian_sparsity
-// row_sparse_Jacobian_sparsity::zero_matrix
-// row_sparse_Jacobian_sparsity::note_new_row
-// row_sparse_Jacobian_sparsity::note_element
+// row_sparse_Jacobian::row_sparse_Jacobian
+// row_sparse_Jacobian::~row_sparse_Jacobian
+// row_sparse_Jacobian::element
+// row_sparse_Jacobian::zero_matrix
+// row_sparse_Jacobian::set_element
+// row_sparse_Jacobian::sum_into_element
+/// row_sparse_Jacobian::find_element
+/// row_sparse_Jacobian::insert_element
+  #ifdef DEBUG_ROW_SPARSE_JACOBIAN
+/// row_sparse_Jacobian::print_data_structure
+  #endif
 #endif
 //
-#ifdef HAVE_ROW_SPARSE_JACOBIAN
-// row_sparse_Jacobian_matrix::row_sparse_Jacobian_matrix
-// row_sparse_Jacobian_matrix::~row_sparse_Jacobian_matrix
-// row_sparse_Jacobian_matrix::zero_matrix
-// row_sparse_Jacobian_matrix::find_element
-// row_sparse_Jacobian_matrix::start_new_row
-// row_sparse_Jacobian_matrix::insert_element
+#ifdef HAVE_ROW_SPARSE_JACOBIAN__ILUCG
+// row_sparse_Jacobian__ILUCG::row_sparse_Jacobian__ILUCG
+// row_sparse_Jacobian__ILUCG::~row_sparse_Jacobian__ILUCG
+// row_sparse_Jacobian__ILUCG::solve_linear_system
 #endif
+//
 
 #include <stdio.h>
 #include <assert.h>
@@ -52,16 +58,102 @@ using jtutil::error_exit;
 //******************************************************************************
 //******************************************************************************
 
+//
+// FIXME:
+//   Cactus's CCTK_FCALL() isn't expanded in .h files (this is a bug),
+//   so we include the contents of "lapack.h" inline here.  This is a
+//   kludge! :( :(
+//#include "ilucg.h"
+//
+
+//***** begin "ilucg.h" contents ******
+/* ilucg.h -- C/C++ prototypes for [sd]ilucg() routines */
+/* $Header$ */
+
+/*
+ * prerequisites:
+ *	"cctk.h"
+ *	"config.h"	// for "integer" = Fortran integer
+ */
+
+#ifdef __cplusplus
+extern "C"
+       {
+#endif
+
+/*
+ * ***** ILUCG *****
+ */
+integer CCTK_FCALL
+  CCTK_FNAME(silucg)(const integer* N,
+		     const integer ia[], const integer ja[], const float a[],
+		     const float b[], float x[],
+		     integer itemp[], float rtemp[],
+		     const float* eps, const integer* iter,
+		     integer* istatus);
+integer CCTK_FCALL
+  CCTK_FNAME(dilucg)(const integer* N,
+		     const integer ia[], const integer ja[], const double a[],
+		     const double b[], double x[],
+		     integer itemp[], double rtemp[],
+		     const double* eps, const integer* iter,
+		     integer* istatus);
+
+/*
+ * ***** wrappers (for returning logical results) *****
+ */
+void CCTK_FCALL
+  CCTK_FNAME(silucg_wrapper)
+		    (const integer* N,
+		     const integer ia[], const integer ja[], const float a[],
+		     const float b[], float x[],
+		     integer itemp[], float rtemp[],
+		     const float* eps, const integer* iter,
+		     integer* istatus, integer* ierror);
+void CCTK_FCALL
+  CCTK_FNAME(dilucg_wrapper)
+		    (const integer* N,
+		     const integer ia[], const integer ja[], const double a[],
+		     const double b[], double x[],
+		     integer itemp[], double rtemp[],
+		     const double* eps, const integer* iter,
+		     integer* istatus, integer* ierror);
+
+#ifdef __cplusplus
+       }	/* extern "C" */
+#endif
+//***** end "ilucg.h" contents ******
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
 #ifdef HAVE_ROW_SPARSE_JACOBIAN
 //
-// This function constructs a  row_sparse_Jacobian_sparsity  object.
+// This function constructs a  row_sparse_Jacobian  object.
 //
-row_sparse_Jacobian_sparsity::row_sparse_Jacobian_sparsity
-	(patch_system& ps,
-	 bool print_msg_flag /* = false */)
-	: Jacobian_sparsity(ps, print_msg_flag),
-	  JJs_at_this_II_(new int[NN()])
+row_sparse_Jacobian::row_sparse_Jacobian(patch_system& ps,
+					 bool print_msg_flag /* = false */)
+	: Jacobian(ps),
+	  N_nonzeros_(0),
+	  N_nonzeros_allocated_(FD_GRID__MOL_AREA * N_rows_),
+				// initial guess, insert_element()
+				// will grow if/when necessary
+	  current_N_rows_(0),
+	  IA_(new integer[N_rows_+1]),
+	  JA_(new integer[N_nonzeros_allocated_]),
+	   A_(new fp     [N_nonzeros_allocated_])
 {
+if (print_msg_flag)
+   then {
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "constructing row_sparse_Jacobian: N_rows_=%d",
+		   N_rows_);
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "   initial N_nonzeros_allocated_=%d",
+		   N_nonzeros_allocated_);
+	}
+
 zero_matrix();
 }
 #endif	// HAVE_ROW_SPARSE_JACOBIAN
@@ -70,11 +162,13 @@ zero_matrix();
 
 #ifdef HAVE_ROW_SPARSE_JACOBIAN
 //
-// This function destroys a  row_sparse_Jacobian_sparsity  object.
+// This function destroys a  row_sparse_Jacobian  object.
 //
-row_sparse_Jacobian_sparsity::~row_sparse_Jacobian_sparsity()
+row_sparse_Jacobian::~row_sparse_Jacobian()
 {
-delete[] JJs_at_this_II_;
+delete[]  A_;
+delete[] JA_;
+delete[] IA_;
 }
 #endif	// HAVE_ROW_SPARSE_JACOBIAN
 
@@ -82,13 +176,13 @@ delete[] JJs_at_this_II_;
 
 #ifdef HAVE_ROW_SPARSE_JACOBIAN
 //
-// This function notes the zeroing of the entire matrix
+// This function gets a matrix element.
 //
-void row_sparse_Jacobian_sparsity::zero_matrix()
+fp row_sparse_Jacobian::element(int II, int JJ)
+	const
 {
-N_nonzeros_ = 0;
-current_II_ = 0;
-N_JJs_at_this_II_ = 0;
+const int fposn = find_element(II,JJ);
+return (fposn > 0) ? fA(fposn) : 0.0;
 }
 #endif	// HAVE_ROW_SPARSE_JACOBIAN
 
@@ -96,12 +190,16 @@ N_JJs_at_this_II_ = 0;
 
 #ifdef HAVE_ROW_SPARSE_JACOBIAN
 //
-// This function notes that we're starting a new row in the matrix traversal.
+// This function zeros a  row_sparse_Jacobian  object.
 //
-void row_sparse_Jacobian_sparsity::note_new_row(int II)
+void row_sparse_Jacobian::zero_matrix()
 {
-current_II_ = II;
-N_JJs_at_this_II_ = 0;
+#ifdef DEBUG_ROW_SPARSE_JACOBIAN
+printf("row_sparse_Jacobian::zero_matrix()\n");
+#endif
+N_nonzeros_ = 0;
+current_N_rows_= 0;
+fIA(1) = 1;
 }
 #endif	// HAVE_ROW_SPARSE_JACOBIAN
 
@@ -109,120 +207,269 @@ N_JJs_at_this_II_ = 0;
 
 #ifdef HAVE_ROW_SPARSE_JACOBIAN
 //
-// This function notes that the (II,JJ) matrix element needs to be stored.
+// This function sets a matrix element to a specified value.
 //
-void row_sparse_Jacobian_sparsity::note_element(int II, int JJ)
+void row_sparse_Jacobian::set_element(int II, int JJ, fp value)
 {
-if (II != current_II_)
-   then {
-	N_nonzeros_ += N_JJs_at_this_II_;	// update totals for old row
-	note_new_row(II);
-	}
-
-// have we already seen this JJ in this row?
-	for (int jindex = 0 ; jindex < N_JJs_at_this_II_ ; ++jindex)
-	{
-	if (JJs_at_this_II_[jindex] == JJ)
-	   then {
-		// yes, we've already seen it
-		// ==> no-op here
-		return;
-		}
-	}
-
-// get to here ==> no, we haven't seen this JJ before in this roW
-//             ==> add it to this row
-if (N_JJs_at_this_II_ >= NN())
-   then error_exit(ERROR_EXIT,
-"***** row_sparse_Jacobian_sparsity():\n"
-"        row buffer overflowed (this should never happen!)\n"
-"        before insertion N_JJs_at_this_II_=%d >= NN()=%d\n"
-"                         N_nonzeros_=%d current_II_=%d\n"
-"        II=%d JJ=%d\n"
-		   ,
-		   N_JJs_at_this_II_, NN(),
-		   N_nonzeros_, current_II_,
-		   II, JJ);					/*NOTREACHED*/
-JJs_at_this_II_[N_JJs_at_this_II_++] = JJ;
+const int fposn = find_element(II,JJ);
+if (fposn > 0)
+   then fA(fposn) = value;
+   else insert_element(II,JJ, value);
 }
 #endif	// HAVE_ROW_SPARSE_JACOBIAN
 
 //******************************************************************************
-//******************************************************************************
-//******************************************************************************
 
+#ifdef HAVE_ROW_SPARSE_JACOBIAN
 //
-// This function constructs a  row_sparse_Jacobian_matrix  object.
+// This function sums a value into a matrix element.
 //
-row_sparse_Jacobian_matrix::row_sparse_Jacobian_matrix
-	(patch_system& ps,
-	 row_sparse_Jacobian_sparsity& JS,
-	 bool print_msg_flag /* = false */)
-	: Jacobian_matrix(ps, JS),
-	  N_nonzeros_(JS.N_nonzeros()),
-	  A_(new fp[N_nonzeros_]),
-	  IA_(new integer[NN()+1]),
-	  JA_(new integer[N_nonzeros_])
+void row_sparse_Jacobian::sum_into_element(int II, int JJ, fp value)
 {
-if (print_msg_flag)
-   then CCTK_VInfo(CCTK_THORNSTRING,
-		  "row_sparse_Jacobian_matrix ctor: N_nonzeros_=%d",
-		  N_nonzeros_);
-
-zero_matrix();
+const int fposn = find_element(II,JJ);
+if (fposn > 0)
+   then fA(fposn) += value;
+   else insert_element(II,JJ, value);
 }
+#endif	// HAVE_ROW_SPARSE_JACOBIAN
 
 //******************************************************************************
 
+#ifdef HAVE_ROW_SPARSE_JACOBIAN
 //
-// This function destroys a  row_sparse_Jacobian_matrix  object.
+// This function searches our data structures for element (II,JJ).
+// If found, it returns the position (1-origin) in A_ and JA_.
+// If not found, it returns 0.
 //
-row_sparse_Jacobian_matrix::~row_sparse_Jacobian_matrix()
+int row_sparse_Jacobian::find_element(int II, int JJ)
+	const
 {
-delete[] JA_;
-delete[] IA_;
-delete[] A_;
+const int fII = fsub(II);
+const int fJJ = fsub(JJ);
+
+if (fII > current_N_rows_)
+   then return 0;				// this row not defined yet
+
+const int fstart = fIA(fII);
+const int fstop  = fIA(fII + 1);
+	for (int fposn = fstart ; fposn < fstop ; ++fposn)
+	{
+	if (fJA(fposn) == fJJ)
+	   then return fposn;				// found
+	}
+
+return 0;						// not found
 }
+#endif	// HAVE_ROW_SPARSE_JACOBIAN
 
 //******************************************************************************
 
+#ifdef HAVE_ROW_SPARSE_JACOBIAN
 //
-// This function zeros a  row_sparse_Jacobian_matrix  object.
+// This function inserts a new element in the matrix, starting a new
+// row and/or growing the arrays if necessary.  It should only be called
+// if JJ isn't already in this row.
 //
-void row_sparse_Jacobian_matrix::zero_matrix()
+int row_sparse_Jacobian::insert_element(int II, int JJ, fp value)
 {
-current_N_nonzeros_ = 0;
-current_II_ = 0;
-JA_[current_II_] = fsub(0);
+const int fII = fsub(II);
+const int fJJ = fsub(JJ);
+
+#ifdef DEBUG_ROW_SPARSE_JACOBIAN
+printf(
+  "row_sparse_Jacobian::insert_element(): fII=%d fJJ=%d\n",
+       fII, fJJ);
+#endif
+
+if (fII != current_N_rows_)
+   then {
+      #ifdef DEBUG_ROW_SPARSE_JACOBIAN
+	printf("   starting new row\n");
+      #endif
+	assert(fII == current_N_rows_+1);
+	assert(current_N_rows_ < N_rows_);
+	++current_N_rows_;
+	fIA(current_N_rows_+1) = fIA(current_N_rows_);
+      #ifdef DEBUG_ROW_SPARSE_JACOBIAN
+	print_data_structure();
+      #endif
+	}
+
+if (fIA(fII+1) > N_nonzeros_allocated_)
+   then {
+      #ifdef DEBUG_ROW_SPARSE_JACOBIAN
+	printf("   growing arrays from N_nonzeros_allocated_=%d",
+	       N_nonzeros_allocated_);
+      #endif
+	N_nonzeros_allocated_ += (N_nonzeros_allocated_ >> 1);
+      #ifdef DEBUG_ROW_SPARSE_JACOBIAN
+	printf(" to %d\n", N_nonzeros_allocated_);
+      #endif
+	integer* const new_JA = new integer[N_nonzeros_allocated_];
+	fp     * const  new_A = new fp     [N_nonzeros_allocated_];
+		for (int posn = 0 ; posn < N_nonzeros_ ; ++posn)
+		{
+		new_JA[posn] = JA_[posn];
+		 new_A[posn] =  A_[posn];
+		}
+	delete[]  A_;
+	delete[] JA_;
+	JA_ = new_JA;
+	 A_ =  new_A;
+	}
+
+++N_nonzeros_;
+const int fposn = fIA(fII+1)++;
+fJA(fposn) = fJJ;
+ fA(fposn) = value;
+#ifdef DEBUG_ROW_SPARSE_JACOBIAN
+printf("   storing at fposn=%d (new N_nonzeros_=%d)\n", fposn, N_nonzeros_);
+#endif
+return fposn;
 }
+#endif	// HAVE_ROW_SPARSE_JACOBIAN
 
 //******************************************************************************
 
+#ifdef HAVE_ROW_SPARSE_JACOBIAN
+#ifdef DEBUG_ROW_SPARSE_JACOBIAN
 //
-// This function searches our data structures for element (II,JJ).
-// If found, it returns the position (0-origin) in A_ and JA_.
-// If not found, it returns -1.
+// This function prints the sparse-matrix data structure.
 //
-int row_sparse_Jacobian_matrix::find_element(int II, int JJ)
+void row_sparse_Jacobian::print_data_structure()
 	const
 {
-const int start = IA(II);
-const int stop = IA(II+1);
-	for (int posn = start ; posn < stop ; ++posn)
+printf("--- begin Jacobian printout\n");
+	for (int fII = 1 ; fII <= current_N_rows_ ; ++fII)
 	{
-	if (JA(posn) == JJ)
-	   then return posn;				// found
+	printf("fII=%d", fII);
+	const int fstart = fIA(fII);
+	const int fstop  = fIA(fII + 1);
+	printf("\tfposn=[%d,%d):", fstart, fstop);
+	printf("\tfJJ =");
+		for (int fposn = fstart ; fposn < fstop ; ++fposn)
+		{
+		printf(" %d", fJA(fposn));
+		}
+	printf("\n");
 	}
+printf("--- end Jacobian printout\n");
+}
+#endif	// DEBUG_ROW_SPARSE_JACOBIAN
+#endif	// HAVE_ROW_SPARSE_JACOBIAN
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+#ifdef HAVE_ROW_SPARSE_JACOBIAN__ILUCG
+//
+// This function constructs a  row_sparse_Jacobian__ILUCG  object.
+//
+row_sparse_Jacobian__ILUCG::row_sparse_Jacobian__ILUCG
+	(patch_system& ps,
+	 bool print_msg_flag /* = false */)
+	: row_sparse_Jacobian(ps, print_msg_flag),
+	  print_msg_flag_(print_msg_flag),
+	  itemp_(NULL),
+	  rtemp_(NULL)
+{ }
+#endif	// HAVE_ROW_SPARSE_JACOBIAN__ILUCG
 
-return -1;						// not found
+//******************************************************************************
+
+#ifdef HAVE_ROW_SPARSE_JACOBIAN__ILUCG
+//
+// This function destroys a  row_sparse_Jacobian__ILUCG  object.
+//
+row_sparse_Jacobian__ILUCG::~row_sparse_Jacobian__ILUCG()
+{
+delete[] rtemp_;
+delete[] itemp_;
 }
+#endif	// HAVE_ROW_SPARSE_JACOBIAN__ILUCG
 
 //******************************************************************************
 
+#ifdef HAVE_ROW_SPARSE_JACOBIAN__ILUCG
 //
-// This function starts a new row in the matrix.
+// This function solves the linear system J.x = rhs, with rhs and x
+// being nominal-grid gridfns, using an ILUCG subroutine originally
+// provided to me in 1985 (!) by Tom Nicol of the UBC Computing Center.
+// It returns -1.0 (no condition number estimate).
 //
-void row_sparse_Jacobian_matrix::start_new_row(int II)
+fp row_sparse_Jacobian__ILUCG::solve_linear_system(int rhs_gfn, int x_gfn)
 {
-// FIXME FIXME
+assert(current_N_rows_ == N_rows_);
+
+//
+// if this is our first call, allocate the scratch arrays
+//
+if (itemp_ == NULL)
+   then {
+	if (print_msg_flag_)
+	   then {
+		CCTK_VInfo(CCTK_THORNSTRING,
+			   "row_sparse_Jacobian__ILUCG::solve_linear_system()");
+		CCTK_VInfo(CCTK_THORNSTRING,
+			   "   N_rows_=%d N_nonzeros_=%d",
+			   N_rows_, N_nonzeros_);
+		CCTK_VInfo(CCTK_THORNSTRING,
+			   "   N_nonzeros_allocated_=%d",
+			   N_nonzeros_allocated_);
+		}
+	itemp_ = new integer[3*N_rows_ + 3*N_nonzeros_ + 2];
+	rtemp_ = new fp     [4*N_rows_ +   N_nonzeros_    ];
+	}
+
+//
+// set up the ILUCG subroutine
+//
+
+// initial guess = all zeros
+fp *x = ps_.gridfn_data(x_gfn);
+	for (int II = 0 ; II < N_rows_ ; ++II)
+	{
+	x[II] = 0.0;
+	}
+
+const integer N = N_rows_;
+const fp *rhs = ps_.gridfn_data(rhs_gfn);
+const fp eps = -1000.0;			// solve to 1000 * machine epsilon
+const integer max_iterations = 0;	// no limit on iterations
+integer istatus, ierror;
+
+// the actual linear solution
+#if   defined(FP_IS_FLOAT)
+  CCTK_FNAME(silucg_wrapper)(&N,
+			     IA_, JA_, A_,
+			     rhs, x,
+			     itemp_, rtemp_,
+			     &eps, &max_iterations,
+			     &istatus, &ierror);
+#elif defined(FP_IS_DOUBLE)
+  CCTK_FNAME(dilucg_wrapper)(&N,
+			     IA_, JA_, A_,
+			     rhs, x,
+			     itemp_, rtemp_,
+			     &eps, &max_iterations,
+			     &istatus, &ierror);
+#else
+  #error "don't know fp datatype!"
+#endif
+
+if (ierror != 0)
+   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"***** row_sparse_Jacobian__ILUCG::solve_linear_system(rhs_gfn=%d, x_gfn=%d):\n"
+"        error return istatus=%d from [sd]ilucg() routine!"
+		   ,
+		   rhs_gfn, x_gfn,
+		   int(istatus));				/*NOTREACHED*/
+if (print_msg_flag_)
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "   solution found in %d CG iterations",
+		   istatus);
+
+return -1.0;			// no condition number estimate available
 }
+#endif	// HAVE_ROW_SPARSE_JACOBIAN__ILUCG