* 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

1 files changed, 33 insertions, 51 deletions

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