* add explicit parameter structure for solve_linear_system()

* change  row_sparse_Jacobian__ILUCG::solve_linear_system()
  to use ILUCG parameters in parameter structure
  for convergence tolerance and iteration limit


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

5 files changed, 70 insertions, 55 deletions

diff --git a/src/elliptic/Jacobian.hh b/src/elliptic/Jacobian.hh
index 0e16182..327cabf 100644
--- a/src/elliptic/Jacobian.hh
+++ b/src/elliptic/Jacobian.hh
@@ -74,6 +74,9 @@
 
 //******************************************************************************
 
+// forward declarations
+class linear_solver_pars;
+
 // ABC to define Jacobian matrix
 class	Jacobian
 	{
@@ -135,6 +138,7 @@ public:
 	//     into the matrix
 	//
 	virtual fp solve_linear_system(int rhs_gfn, int x_gfn,
+				       const struct linear_solver_pars& pars,
 				       bool print_msg_flag) = 0;
 
 
@@ -162,6 +166,24 @@ protected:
 	int N_rows_;
 	};
 
+//**************************************
+
+//
+// this class defines parameters for  solve_linear_system()
+// for all of our derived classes
+//
+struct	linear_solver_pars
+	{
+	struct	LAPACK_pars
+		{
+		} LAPACK_pars;
+	struct	ILUCG_pars
+		{
+		fp  error_tolerance;
+		int max_CG_iterations;
+		} ILUCG_pars;
+	};
+
 //******************************************************************************
 
 enum	Jacobian_store_solve_method
diff --git a/src/elliptic/dense_Jacobian.cc b/src/elliptic/dense_Jacobian.cc
index 3e4e582..4fd9e66 100644
--- a/src/elliptic/dense_Jacobian.cc
+++ b/src/elliptic/dense_Jacobian.cc
@@ -165,8 +165,9 @@ void dense_Jacobian::zero_matrix()
 //
 // This function constructs a  dense_Jacobian__LAPACK  object.
 //
-dense_Jacobian__LAPACK::dense_Jacobian__LAPACK(patch_system& ps,
-			       bool print_msg_flag /* = false */)
+dense_Jacobian__LAPACK::dense_Jacobian__LAPACK
+	(patch_system& ps,
+	 bool print_msg_flag /* = false */)
 	: dense_Jacobian(ps, print_msg_flag),
 	  pivot_(new integer[  N_rows_]),
 	  iwork_(new integer[  N_rows_]),
@@ -201,8 +202,10 @@ delete[] pivot_;
 // It returns the estimated infinity-norm condition number of the linear
 // system, or 0.0 if the matrix is numerically singular
 //
-fp dense_Jacobian__LAPACK::solve_linear_system(int rhs_gfn, int x_gfn,
-					       bool print_msg_flag)
+fp dense_Jacobian__LAPACK::solve_linear_system
+	(int rhs_gfn, int x_gfn,
+	 const struct linear_solver_pars& pars,
+	 bool print_msg_flag)
 {
 const fp *rhs = ps_.gridfn_data(rhs_gfn);
 fp       *x   = ps_.gridfn_data(x_gfn);
diff --git a/src/elliptic/dense_Jacobian.hh b/src/elliptic/dense_Jacobian.hh
index 883f2ea..6946993 100644
--- a/src/elliptic/dense_Jacobian.hh
+++ b/src/elliptic/dense_Jacobian.hh
@@ -91,6 +91,7 @@ public:
 	//	       0.0 if matrix is numerically singular,
 	//             -1.0 if condition number is unknown
 	fp solve_linear_system(int rhs_gfn, int x_gfn,
+			       const struct linear_solver_pars& pars,
 			       bool print_msg_flag);
 
 	// constructor, destructor
@@ -108,6 +109,12 @@ private:
 	integer *iwork_;	// size N_rows_
 	fp *rwork_;		// size 4*N_rows_
 	};
+
+//**************************************
+
+struct	Jacobian_pars__LAPACK
+	: public Jacobian_pars;
+
 #endif	// HAVE_DENSE_JACOBIAN__LAPACK
 
 //******************************************************************************
diff --git a/src/elliptic/row_sparse_Jacobian.cc b/src/elliptic/row_sparse_Jacobian.cc
index 6c292cb..8dc6ff4 100644
--- a/src/elliptic/row_sparse_Jacobian.cc
+++ b/src/elliptic/row_sparse_Jacobian.cc
@@ -61,9 +61,9 @@ 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"
+//   so we include the contents of "../sparse-matrix/ilucg/ilucg.h"
+//   inline here.  This is a doubleplusungood kludge! :( :(
+//#include "../sparse-matrix/ilucg/ilucg.h"
 //
 
 //***** begin "ilucg.h" contents ******
@@ -84,40 +84,20 @@ extern "C"
 /*
  * ***** 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);
+  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* 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);
+  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, integer* IERROR);
 
 #ifdef __cplusplus
        }	/* extern "C" */
@@ -397,8 +377,10 @@ delete[] itemp_;
 // provided to me in 1985 (!) by Tom Nicol of the UBC Computing Center.
 // It returns -1.0 (no condition number estimate).
 //
-fp row_sparse_Jacobian__ILUCG::solve_linear_system(int rhs_gfn, int x_gfn,
-						   bool print_msg_flag)
+fp row_sparse_Jacobian__ILUCG::solve_linear_system
+	(int rhs_gfn, int x_gfn,
+	 const struct linear_solver_pars& pars,
+	 bool print_msg_flag)
 {
 assert(current_N_rows_ == N_rows_);
 
@@ -435,25 +417,25 @@ fp *x = ps_.gridfn_data(x_gfn);
 
 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
+const fp      eps            = pars.ILUCG_pars.error_tolerance;
+const integer max_iterations = pars.ILUCG_pars.max_CG_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);
+  CCTK_FNAME(silucg)(&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);
+  CCTK_FNAME(dilucg)(&N,
+		     IA_, JA_, A_,
+		     rhs, x,
+		     itemp_, rtemp_,
+		     &eps, &max_iterations,
+		     &istatus, &ierror);
 #else
   #error "don't know fp datatype!"
 #endif
diff --git a/src/elliptic/row_sparse_Jacobian.hh b/src/elliptic/row_sparse_Jacobian.hh
index 1c11878..eae865e 100644
--- a/src/elliptic/row_sparse_Jacobian.hh
+++ b/src/elliptic/row_sparse_Jacobian.hh
@@ -189,6 +189,7 @@ public:
 	//	       0.0 if matrix is numerically singular,
 	//             -1.0 if condition number is unknown
 	fp solve_linear_system(int rhs_gfn, int x_gfn,
+			       const struct linear_solver_pars& pars,
 			       bool print_msg_flag);
 
 	// constructor, destructor