From 8b0a74a327fbecb12f3d32ecff1fa28a6a6c2a15 Mon Sep 17 00:00:00 2001 From: jthorn Date: Wed, 22 Jan 2003 12:27:46 +0000 Subject: * 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 --- src/elliptic/row_sparse_Jacobian.cc | 455 +++++++++++++++++++++++++++--------- 1 file changed, 351 insertions(+), 104 deletions(-) (limited to 'src/elliptic/row_sparse_Jacobian.cc') 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 @@ -1,23 +1,29 @@ // Jacobian.cc -- data structures for the Jacobian matrix // $Header$ +// +// <<>> // #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 #include @@ -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 -- cgit v1.2.3