From b1fb8f18cc2231bde4ea6d5d71348bd7894af370 Mon Sep 17 00:00:00 2001 From: jthorn Date: Mon, 2 Jun 2003 15:11:57 +0000 Subject: * 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 --- src/elliptic/Jacobian.hh | 22 ++++++++++ src/elliptic/dense_Jacobian.cc | 11 +++-- src/elliptic/dense_Jacobian.hh | 7 ++++ src/elliptic/row_sparse_Jacobian.cc | 84 +++++++++++++++---------------------- src/elliptic/row_sparse_Jacobian.hh | 1 + 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 -- cgit v1.2.3