diff options
author | jthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5> | 2003-01-20 14:01:12 +0000 |
---|---|---|
committer | jthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5> | 2003-01-20 14:01:12 +0000 |
commit | 4fc46cfe1961ba2b7c54862f561d48f87e10a116 (patch) | |
tree | 1c10d76ff23e249765d9ba943fb4a8e1a48e77f8 | |
parent | f595463096694fd5901bf2503d20d0e30269653d (diff) |
start adding support for sparse matrix Jacobians (doesn't work yet)
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@926 f88db872-0e4f-0410-b76b-b9085cfa78c5
-rw-r--r-- | src/driver/Newton.cc | 2 | ||||
-rw-r--r-- | src/driver/driver.hh | 13 | ||||
-rw-r--r-- | src/driver/find_horizons.cc | 42 | ||||
-rw-r--r-- | src/driver/io.cc | 41 | ||||
-rw-r--r-- | src/driver/setup.cc | 7 | ||||
-rw-r--r-- | src/elliptic/Jacobian.cc | 300 | ||||
-rw-r--r-- | src/elliptic/Jacobian.hh | 228 | ||||
-rw-r--r-- | src/elliptic/make.code.defn | 4 | ||||
-rw-r--r-- | src/gr/expansion_Jacobian.cc | 40 | ||||
-rw-r--r-- | src/gr/gr.hh | 8 | ||||
-rw-r--r-- | src/gr/misc.cc | 47 | ||||
-rw-r--r-- | src/include/config.h | 6 |
12 files changed, 337 insertions, 401 deletions
diff --git a/src/driver/Newton.cc b/src/driver/Newton.cc index a93e101..622bece 100644 --- a/src/driver/Newton.cc +++ b/src/driver/Newton.cc @@ -60,7 +60,7 @@ using jtutil::error_exit; // otherwise it's false. // bool Newton_solve(patch_system& ps, - Jacobian& Jac, + Jacobian_matrix& Jac, const struct cactus_grid_info& cgi, const struct geometry_info& gi, const struct Jacobian_info& Jacobian_info, diff --git a/src/driver/driver.hh b/src/driver/driver.hh index 119bafe..9747002 100644 --- a/src/driver/driver.hh +++ b/src/driver/driver.hh @@ -181,7 +181,7 @@ struct BH_diagnostics struct AH_info { patch_system* ps_ptr; - Jacobian* Jac_ptr; + Jacobian_matrix* Jac_ptr; struct initial_guess_info initial_guess_info; @@ -245,7 +245,7 @@ enum initial_guess_method // Newton.cc // returns true for success, false for failure to converge bool Newton_solve(patch_system& ps, - Jacobian& Jac, + Jacobian_matrix& Jac, const struct cactus_grid_info& cgi, const struct geometry_info& gi, const struct Jacobian_info& Jacobian_info, @@ -262,10 +262,11 @@ void input_gridfn(patch_system& ps, int unknown_gfn, void output_gridfn(patch_system& ps, int unknown_gfn, const struct IO_info& IO_info, const char base_file_name[], int hn, bool print_msg_flag, int AHF_iteration = 0); -void print_Jacobians(const patch_system& ps, - const Jacobian* Jac_NP, const Jacobian* Jac_SD_FDdr, - const struct IO_info& IO_info, const char base_file_name[], - int hn, bool print_msg_flag, int AHF_iteration = 0); +void output_Jacobians(const patch_system& ps, + const Jacobian_matrix* Jac_NP_ptr, + const Jacobian_matrix* Jac_SD_FDdr_ptr, + const struct IO_info& IO_info, const char base_file_name[], + int hn, bool print_msg_flag, int AHF_iteration = 0); FILE* setup_BH_diagnostics_output_file(const struct IO_info& IO_info, int hn, int N_horizons); void output_BH_diagnostics_fn(const struct BH_diagnostics& BH_diagnostics, diff --git a/src/driver/find_horizons.cc b/src/driver/find_horizons.cc index e721041..45710dd 100644 --- a/src/driver/find_horizons.cc +++ b/src/driver/find_horizons.cc @@ -80,7 +80,7 @@ bool do_test_expansion_Jacobian bool test_all_Jacobian_methods, struct Jacobian_info& Jac_info, struct cactus_grid_info& cgi, struct geometry_info& gi, - patch_system& ps, Jacobian* Jac_ptr, + patch_system& ps, Jacobian_matrix& Jac, bool active_flag, int hn, int N_horizons); bool do_find_horizon (const struct verbose_info& verbose_info, int timer_handle, @@ -88,7 +88,7 @@ bool do_find_horizon const struct solver_info& solver_info, bool initial_find_flag, const struct Jacobian_info& Jac_info, struct cactus_grid_info& cgi, struct geometry_info& gi, - patch_system& ps, Jacobian* Jac_ptr, + patch_system& ps, Jacobian_matrix& Jac, bool active_flag, int hn, int N_horizons); void compute_BH_diagnostics @@ -179,6 +179,16 @@ if (state.gi.geometry_method == geometry__local_interp_from_Cactus_grid) .print_algorithm_highlights); } + // if we're going to need a Jacobian + // and it doesn't exist, create it. + if ( (state.method != method__evaluate_expansion) + && (AH_info.Jac_ptr == NULL) ) + then AH_info.Jac_ptr + = create_Jacobian_matrix(ps, + state.cgi, state.gi, + state.Jac_info, + verbose_info.print_algorithm_details); + AH_info.AH_found = false; switch (state.method) { @@ -195,7 +205,7 @@ if (state.gi.geometry_method == geometry__local_interp_from_Cactus_grid) IO_info, test_all_Jacobian_methods, state.Jac_info, state.cgi, state.gi, - ps, AH_info.Jac_ptr, + ps, *AH_info.Jac_ptr, active_flag, hn, state.N_horizons); break; @@ -205,7 +215,7 @@ if (state.gi.geometry_method == geometry__local_interp_from_Cactus_grid) IO_info, output_h, output_Theta, solver_info, initial_find_flag, state.Jac_info, state.cgi, state.gi, - ps, AH_info.Jac_ptr, + ps, *AH_info.Jac_ptr, active_flag, hn, state.N_horizons); // store found flag in Cactus array variable @@ -458,11 +468,11 @@ bool do_test_expansion_Jacobian bool test_all_Jacobian_methods, struct Jacobian_info& Jac_info, struct cactus_grid_info& cgi, struct geometry_info& gi, - patch_system& ps, Jacobian* Jac_ptr, + patch_system& ps, Jacobian_matrix& Jac, bool active_flag, int hn, int N_horizons) { // numerical perturbation -Jacobian* Jac_NP_ptr = Jac_ptr; +Jacobian_matrix* Jac_NP_ptr = & Jac; if (! expansion(ps, cgi, gi, true)) then return false; // *** ERROR RETURN *** Jac_info.Jacobian_method = Jacobian_method__numerical_perturb; @@ -471,11 +481,15 @@ if (! expansion_Jacobian(ps, *Jac_NP_ptr, true)) // print msgs then return false; // *** ERROR RETURN *** -Jacobian* Jac_SD_FDdr_ptr = NULL; +Jacobian_matrix* Jac_SD_FDdr_ptr = NULL; if (test_all_Jacobian_methods) then { // symbolic differentiation with finite diff d/dr - Jac_SD_FDdr_ptr = & new_Jacobian(ps, Jac_info.Jacobian_storage_method); + Jac_SD_FDdr_ptr + = create_Jacobian_matrix(ps, + cgi, gi, + Jac_info, + verbose_info.print_algorithm_details); if (! expansion(ps, cgi, gi, true)) then return false; // *** ERROR RETURN *** Jac_info.Jacobian_method = Jacobian_method__symbolic_diff_with_FD_dr; @@ -486,10 +500,10 @@ if (test_all_Jacobian_methods) } if (active_flag) - then print_Jacobians(ps, - Jac_NP_ptr, Jac_SD_FDdr_ptr, - IO_info, IO_info.Jacobian_base_file_name, - hn, true); + then output_Jacobians(ps, + Jac_NP_ptr, Jac_SD_FDdr_ptr, + IO_info, IO_info.Jacobian_base_file_name, + hn, true); return true; // *** NORMAL RETURN *** } @@ -534,11 +548,9 @@ bool do_find_horizon const struct solver_info& solver_info, bool initial_find_flag, const struct Jacobian_info& Jac_info, struct cactus_grid_info& cgi, struct geometry_info& gi, - patch_system& ps, Jacobian* Jac_ptr, + patch_system& ps, Jacobian_matrix& Jac, bool active_flag, int hn, int N_horizons) { -Jacobian& Jac = *Jac_ptr; - if (verbose_info.print_algorithm_highlights) then CCTK_VInfo(CCTK_THORNSTRING, "searching for horizon %d/%d", diff --git a/src/driver/io.cc b/src/driver/io.cc index cc6be75..23c6fa2 100644 --- a/src/driver/io.cc +++ b/src/driver/io.cc @@ -2,9 +2,9 @@ // $Header$ // // <<<prototypes for functions local to this file>>> -// input_gridfn - read an angular grid function from a data file -// output_gridfn - write an angular grid function to a data file -// output_Jacobians - write a Jacobian matrix or matrices to a data file +// input_gridfn - read an angular grid function from an input file +// output_gridfn - write an angular grid function to an output file +// output_Jacobians - write a Jacobian matrix or matrices to an output file // decode_horizon_file_format - decode the horizon_file_format parameter // // setup_BH_diagnostics_output_file - create/initialize BH-diagnostics file @@ -186,15 +186,16 @@ default: // Arguments: // A null Jacobian pointer means to skip that Jacobian. // -void print_Jacobians(const patch_system& ps, - const Jacobian* Jac_NP, const Jacobian* Jac_SD_FDdr, - const struct IO_info& IO_info, const char base_file_name[], - int hn, bool print_msg_flag, int AHF_iteration /* = 0 */) +void output_Jacobians(const patch_system& ps, + const Jacobian_matrix* Jac_NP_ptr, + const Jacobian_matrix* Jac_SD_FDdr_ptr, + const struct IO_info& IO_info, const char base_file_name[], + int hn, bool print_msg_flag, int AHF_iteration /* = 0 */) { -if (Jac_NP == NULL) +if (Jac_NP_ptr == NULL) then { CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, - "print_Jacobinas(): Jac_NP == NULL is not (yet) supported!"); + "output_Jacobians(): Jac_NP_ptr == NULL is not (yet) supported!"); return; // *** ERROR RETURN *** } @@ -203,14 +204,14 @@ const char* file_name = io_file_name(IO_info, base_file_name, if (print_msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, " writing %s to \"%s\"", - ((Jac_SD_FDdr == NULL) ? "NP Jacobian" - : "NP and SD_FDdr Jacobians"), + ((Jac_SD_FDdr_ptr == NULL) ? "NP Jacobian" + : "NP and SD_FDdr Jacobians"), file_name); FILE *fileptr = fopen(file_name, "w"); if (fileptr == NULL) then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, - "print_Jacobians(): can't open output file \"%s\"!", + "output_Jacobians(): can't open output file \"%s\"!", file_name); /*NOTREACHED*/ fprintf(fileptr, "# column 1 = x patch name\n"); @@ -223,10 +224,10 @@ fprintf(fileptr, "# column 7 = y patch number\n"); fprintf(fileptr, "# column 8 = y irho\n"); fprintf(fileptr, "# column 9 = y isigma\n"); fprintf(fileptr, "# column 10 = y JJ\n"); -fprintf(fileptr, "# column 11 = Jac_NP(II,JJ)\n"); -if (Jac_SD_FDdr != NULL) +fprintf(fileptr, "# column 11 = Jac_NP.element(II,JJ)\n"); +if (Jac_SD_FDdr_ptr != NULL) then { - fprintf(fileptr, "# column 12 = Jac_SD_FDdr(II,JJ)\n"); + fprintf(fileptr, "# column 12 = Jac_SD_FDdr.element(II,JJ)\n"); fprintf(fileptr, "# column 13 = abs error in Jac_SD_FDdr\n"); fprintf(fileptr, "# column 14 = rel error in Jac_SD_FDdr\n"); } @@ -257,11 +258,13 @@ if (Jac_SD_FDdr != NULL) { const int JJ = ps.gpn_of_patch_irho_isigma(yp, y_irho,y_isigma); - if (! Jac_NP->is_explicitly_stored(II,JJ)) + if (! Jac_NP_ptr->is_explicitly_stored(II,JJ)) then continue; // skip sparse points - const fp NP = (*Jac_NP)(II,JJ); - const fp SD_FDdr = (Jac_SD_FDdr == NULL) ? 0.0 : (*Jac_SD_FDdr)(II,JJ); + const fp NP = Jac_NP_ptr->element(II,JJ); + const fp SD_FDdr = (Jac_SD_FDdr_ptr == NULL) + ? 0.0 + : Jac_SD_FDdr_ptr->element(II,JJ); if ((NP == 0.0) && (SD_FDdr == 0.0)) then continue; // skip zero values (if == ) @@ -272,7 +275,7 @@ if (Jac_SD_FDdr != NULL) yp.name(), ypn, y_irho, y_isigma, JJ, double(NP)); - if (Jac_SD_FDdr != NULL) + if (Jac_SD_FDdr_ptr != NULL) then { const fp abs_NP = jtutil::abs(NP ); const fp abs_SD_FDdr = jtutil::abs(SD_FDdr); diff --git a/src/driver/setup.cc b/src/driver/setup.cc index afebbdd..ad4e971 100644 --- a/src/driver/setup.cc +++ b/src/driver/setup.cc @@ -311,12 +311,7 @@ state.AH_info_array = new AH_info[N_horizons+1]; ps.set_gridfn_to_constant(1.0, gfns::gfn__one); // create the Jacobian matrix - if (verbose_info.print_algorithm_details) - then CCTK_VInfo(CCTK_THORNSTRING, - " constructing Jacobian matrix"); - AH_info.Jac_ptr = (state.method == method__evaluate_expansion) - ? NULL // no Jacobian used for this case - : &new_Jacobian(ps, Jac_info.Jacobian_storage_method); + AH_info.Jac_ptr = NULL; if (verbose_info.print_algorithm_details) then CCTK_VInfo(CCTK_THORNSTRING, diff --git a/src/elliptic/Jacobian.cc b/src/elliptic/Jacobian.cc index 4ac18c3..048c6ac 100644 --- a/src/elliptic/Jacobian.cc +++ b/src/elliptic/Jacobian.cc @@ -1,21 +1,12 @@ -// Jacobian.cc -- data structures for the Jacobian matrix +// Jacobian.cc -- generic routines for the Jacobian matrix // $Header$ // // <<<literal contents of "lapack.hh">>> // // decode_Jacobian_type -// new_Jacobian -// -// Jacobian::Jacobian -// -#ifdef HAVE_DENSE_JACOBIAN -// dense_Jacobian::dense_Jacobian -// dense_Jacobian::~dense_Jacobian -// dense_Jacobian::zero_matrix -// dense_Jacobian::zero_row -// dense_Jacobian::solve_linear_system -#endif +// new_Jacobian_sparsity +// new_Jacobian_matrix // #include <stdio.h> @@ -45,76 +36,8 @@ using jtutil::error_exit; #include "../patch/patch_system.hh" #include "Jacobian.hh" - -//****************************************************************************** -//****************************************************************************** -//****************************************************************************** - -// -// 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. -//#include "lapack.h" -// - -//***** begin "lapack.h" contents ****** -/* lapack.h -- C/C++ prototypes for (some) BLAS+LAPACK+wrapper routines */ -/* $Header$ */ - -/* - * prerequisites: - * "cctk.h" - * "config.hh" - */ - -#ifdef __cplusplus -extern "C" { -#endif - -/* - * ***** BLAS ***** - */ -integer CCTK_FCALL - CCTK_FNAME(isamax)(const integer* N, const float SX[], const integer* incx); -integer CCTK_FCALL - CCTK_FNAME(idamax)(const integer* N, const double DX[], const integer* incx); - -/* - * ***** LAPACK ***** - */ -void CCTK_FCALL - CCTK_FNAME(sgesv)(const integer* N, const integer* NRHS, - float A[], const int* LDA, - integer IPIV[], - float B[], const integer* LDB, integer* info); -void CCTK_FCALL - CCTK_FNAME(dgesv)(const integer* N, const integer* NRHS, - double A[], const int* LDA, - integer IPIV[], - double B[], const integer* LDB, integer* info); - -/* - * ***** wrappers (for passing character-string args) ***** - */ -/* norm_int = 0 for infinity-norm, 1 for 1-norm */ -void CCTK_FCALL - CCTK_FNAME(sgecon_wrapper)(const integer* norm_int, - const integer* N, - float A[], const integer* LDA, - const float* anorm, float* rcond, - float WORK[], integer IWORK[], - integer* info); -void CCTK_FCALL - CCTK_FNAME(dgecon_wrapper)(const integer* norm_int, - const integer* N, - double A[], const integer* LDA, - const double* anorm, double* rcond, - double WORK[], integer IWORK[], - integer* info); - -#ifdef __cplusplus - } /* extern "C" */ -#endif -//***** end "lapack.h" contents ******** +#include "dense_Jacobian.hh" +#include "row_sparse_Jacobian.hh" //****************************************************************************** //****************************************************************************** @@ -131,200 +54,101 @@ if (STRING_EQUAL(Jacobian_type_string, "dense matrix")) then { #ifdef HAVE_DENSE_JACOBIAN return Jacobian_type__dense_matrix; - #else - error_exit(ERROR_EXIT, -"\n" -" decode_Jacobian_type():\n" -" Jacobian type \"dense matrix\" is not configured in this binary;\n" -" see \"src/include/config.hh\" for details on what Jacobian types\n" -" are configured and how to change this\n"); /*NOTREACHED*/ #endif } + +else if (STRING_EQUAL(Jacobian_type_string, "row-oriented sparse matrix")) + then { + #ifdef HAVE_ROW_SPARSE_JACOBIAN + return Jacobian_type__row_sparse_matrix; + #endif + } + else error_exit(ERROR_EXIT, "decode_Jacobian_type(): unknown Jacobian_type_string=\"%s\"!\n", Jacobian_type_string); /*NOTREACHED*/ + +// fall through to here ==> we recognize the matrix type, +// but it's not configured +error_exit(ERROR_EXIT, +"\n" +" decode_Jacobian_type():\n" +" Jacobian type \"%s\" is not configured in this binary;\n" +" see \"src/include/config.hh\" for details on what Jacobian types\n" +" are configured and how to change this\n" + , + Jacobian_type_string); /*NOTREACHED*/ } //****************************************************************************** // -// This function is an "object factory" for Jacobians: it constructs -// and returns a new-allocated Jacobian object of a specified type. +// This function is an "object factory" for Jacobian_sparsity objects: +// it constructs and returns a pointer to a new-allocated Jacobian_sparsity +// object of the specified derived type. // // FIXME: the patch system shouldn't really have to be non-const, but // the Jacobian constructors all require this to allow the // linear solvers to directly update gridfns // -Jacobian& new_Jacobian(patch_system& ps, enum Jacobian_type Jac_type) +Jacobian_sparsity* new_Jacobian_sparsity(enum Jacobian_type Jac_type, + patch_system& ps, + bool print_msg_flag /* = false */) { switch (Jac_type) { #ifdef HAVE_DENSE_JACOBIAN case Jacobian_type__dense_matrix: - return *new dense_Jacobian(ps); + return new dense_Jacobian_sparsity(ps, print_msg_flag); +#endif + +#ifdef HAVE_ROW_SPARSE_JACOBIAN + case Jacobian_type__row_sparse_matrix: + return new row_sparse_Jacobian_sparsity(ps, print_msg_flag); #endif + default: error_exit(ERROR_EXIT, -"new_Jacobian(): unknown Jacobian_type=(int)%d!\n", +"new_Jacobian_sparsity(): unknown Jacobian_type=(int)%d!\n", Jac_type); /*NOTREACHED*/ } } //****************************************************************************** -//****************************************************************************** -//****************************************************************************** // -// This function constructs a Jacobian object. +// This function is an "object factory" for Jacobian_matrix objects: +// it constructs and returns a pointer to a new-allocated Jacobian_matrix +// object of the specified derived type. // -Jacobian::Jacobian(patch_system& ps) - : ps_(ps), - NN_(ps.N_grid_points()) -{ } - -//****************************************************************************** -//****************************************************************************** -//****************************************************************************** - -#ifdef HAVE_DENSE_JACOBIAN -// -// This function constructs a dense_Jacobian object. -// -dense_Jacobian::dense_Jacobian(patch_system& ps) - : Jacobian(ps), - matrix_(0,NN()-1, 0,NN()-1), - pivot_(new integer[NN()]), - iwork_(new integer[NN()]), - rwork_(new fp[4*NN()]) // no comma -{ } -#endif - -//****************************************************************************** - -#ifdef HAVE_DENSE_JACOBIAN -// -// THis function destroys a dense_Jacobian object. -// -dense_Jacobian::~dense_Jacobian() -{ -delete[] iwork_; -delete[] rwork_; -delete[] pivot_; -} -#endif - -//****************************************************************************** - -#ifdef HAVE_DENSE_JACOBIAN -// -// This function zeros a dense_Jacobian object. +// FIXME: the patch system shouldn't really have to be non-const, but +// the Jacobian constructors all require this to allow the +// linear solvers to directly update gridfns // -void dense_Jacobian::zero_matrix() +Jacobian_matrix* new_Jacobian_matrix(enum Jacobian_type Jac_type, + patch_system& ps, + Jacobian_sparsity& JS, + bool print_msg_flag /* = false */) { - for (int JJ = 0 ; JJ < NN() ; ++JJ) +switch (Jac_type) { - for (int II = 0 ; II < NN() ; ++II) - { - operator()(II,JJ) = 0.0; - } - } -} -#endif - -//****************************************************************************** - #ifdef HAVE_DENSE_JACOBIAN -// -// This function zeros a single row of a dense_Jacobian object. -// -void dense_Jacobian::zero_row(int II) -{ - for (int JJ = 0 ; JJ < NN() ; ++JJ) - { - operator()(II,JJ) = 0.0; - } -} + case Jacobian_type__dense_matrix: + return new dense_Jacobian_matrix + (ps, static_cast<dense_Jacobian_sparsity&>(JS), + print_msg_flag); #endif -//****************************************************************************** - -#ifdef HAVE_DENSE_JACOBIAN -// -// This function solves the linear system J.x = rhs, with rhs and x -// being nominal-grid gridfns, using LAPACK LU-decomposition routines. -// It returns the estimated infinity-norm condition number of the linear -// system, or 0.0 if the matrix is numerically singular -// -fp dense_Jacobian::solve_linear_system(int rhs_gfn, int x_gfn) -{ -const fp *rhs = my_patch_system().gridfn_data(rhs_gfn); -fp *x = my_patch_system().gridfn_data(x_gfn); -fp *A = matrix_.data_array(); - -const integer infinity_norm_flag = 0; -const integer stride = 1; -const integer NRHS = 1; -const integer N = NN(); -const integer N2 = NN() * NN(); - -// compute the infinity-norm of the matrix A -// ... max_posn = 1-origin index of A[] element with largest absolute value -#if defined(FP_IS_FLOAT) - const integer max_posn = CCTK_FNAME(isamax)(&N2, A, &stride); -#elif defined(FP_IS_DOUBLE) - const integer max_posn = CCTK_FNAME(idamax)(&N2, A, &stride); -#else - #error "don't know fp datatype!" +#ifdef HAVE_ROW_SPARSE_JACOBIAN + case Jacobian_type__row_sparse_matrix: + return new row_sparse_Jacobian_matrix + (ps, static_cast<row_sparse_Jacobian_sparsity&>(JS), + print_msg_flag); #endif -const fp A_infnorm = jtutil::abs(A[max_posn-1]); -// LU decompose and solve the linear system -// -// ... [sd]gesv() use an "in out" design, where the same argument -// is used for both rhs and x ==> we must first copy rhs to x -// - for (int II = 0 ; II < NN() ; ++II) - { - x[II] = rhs[II]; + default: + error_exit(ERROR_EXIT, +"new_Jacobian_matrix(): unknown Jacobian_type=(int)%d!\n", + Jac_type); /*NOTREACHED*/ } -integer info; -#if defined(FP_IS_FLOAT) - CCTK_FNAME(sgesv)(&N, &NRHS, A, &N, pivot_, x, &N, &info); -#elif defined(FP_IS_DOUBLE) - CCTK_FNAME(dgesv)(&N, &NRHS, A, &N, pivot_, x, &N, &info); -#else - #error "don't know fp datatype!" -#endif - -if (info < 0) - then error_exit(ERROR_EXIT, -"\n" -"***** dense_Jacobian::solve_linear_system(rhs_gfn=%d, x_gfn=%d):\n" -" error return (bad argument) info=%d from [sd]gesv() LAPACK routine!\n" -, - rhs_gfn, x_gfn, - int(info)); /*NOTREACHED*/ - -if (info > 0) - then return 0.0; // *** ERROR RETURN *** - // *** (singular matrix) - -// estimate infinity-norm condition number -fp rcond; -#if defined(FP_IS_FLOAT) - CCTK_FNAME(sgecon_wrapper)(&infinity_norm_flag, - &N, A, &N, &A_infnorm, &rcond, - rwork_, iwork_, &info); -#elif defined(FP_IS_DOUBLE) - CCTK_FNAME(dgecon_wrapper)(&infinity_norm_flag, - &N, A, &N, &A_infnorm, &rcond, - rwork_, iwork_, &info); -#else - #error "don't know fp datatype!" -#endif -if (rcond == 0.0) - then return 0.0; // *** ERROR RETURN *** - // *** (singular matrix) -return 1.0/rcond; } -#endif diff --git a/src/elliptic/Jacobian.hh b/src/elliptic/Jacobian.hh index 3faf2d1..d4b90ab 100644 --- a/src/elliptic/Jacobian.hh +++ b/src/elliptic/Jacobian.hh @@ -1,11 +1,10 @@ -// Jacobian.hh -- data structures for the Jacobian matrix +// Jacobian.hh -- generic data structures for the Jacobian matrix // $Header$ // -// Jacobian -- abstract base class to describe a Jacobian matrix -#ifdef HAVE_DENSE_JACOBIAN -// dense_Jacobian -- Jacobian stored as a dense matrix -#endif +// Jacobian -- ABC to describe Jacobian matrix +// Jacobian_sparsity -- ABC to accumulate Jacobian sparsity info +// Jacobian_matrix -- ABC to store/use Jacobian matrix // // decode_Jacobian_type - decode string into internal enum // new_Jacobian - factory method @@ -13,29 +12,53 @@ // // prerequisites: -// "stdc.h" -// "config.hh" -// "../jtutil/util.hh" -// "../jtutil/array.hh" -// "../jtutil/cpm_map.hh" -// "../jtutil/linear_map.hh" -// "coords.hh" -// "grid.hh" -// "fd_grid.hh" -// "patch.hh" -// "patch_edge.hh" -// "patch_interp.hh" -// "ghost_zone.hh" -// "patch_system.hh" +// "../patch/patch_system.hh" // //****************************************************************************** +//****************************************************************************** +//****************************************************************************** // -// A Jacobian object stores the (a) Jacobian matrix for a patch system. -// Jacobian is an abstract base class, from which we derive a separate -// concrete class for each type of sparsity pattern (more precisely for -// each type of storage scheme) of the Jacobian matrix. +// These classes are used to store and manipulate Jacobian matrices +// for a patch system. +// +// The inheritance diagram for the Jacobian classes looks like this: +// Jacobian +// Jacobian_sparsity +// dense_Jacobian_sparsity #ifdef HAVE_DENSE_JACOBIAN +// row_sparse_Jacobian_sparsity #ifdef HAVE_ROW_SPARSE_JACOBIAN +// Jacobian_matrix +// dense_Jacobian_matrix #ifdef HAVE_DENSE_JACOBIAN +// row_sparse_Jacobian_matrix #ifdef HAVE_ROW_SPARSE_JACOBIAN +// where the most-derived classes are all conditional on the noted #ifdef +// symbols. +// +// The Jacobian_sparsity objects are used to accumulate information +// on the Jacobian's sparsity pattern; this is then used in constructing +// the corresponding Jacobian_matrix object (which actually stores +// the Jacobian, and contains member functions to solve elliptic systems +// using the Jacobian as the LHS matrix). +// +// We assume a "traversal" interface for computing Jacobians, defined +// by the Jacobian API. Because this API is shared by (all) the +// Jacobian_sparsity and Jacobian_matrix classes, the same traversal +// routines (prototyped to update Jacobian objects) can be used both +// to record sparsity patterns (ignoring the actual floating-point values) +// and to store the Jacobian matrix. +// +// Thus the typical code to construct and use a Jacobian matrix looks +// like this: +// // prototype +// void trverse_Jacobian(const patch_system& ps, Jacobian& Jac); +// +// Jacobian_sparsity& JS = new_Jacobian_sparsity(Jac_type, ps); +// traverse_Jacobian(ps, JS); +// Jacobian_matrix& Jac = new_Jacobian_matrix(Jac_type, ps, JS); +// traverse_Jacobian(ps, Jac); +// Jac.solve_linear_system(...); +// + // // A row/column index of the Jacobian (denoted II/JJ) is a 0-origin grid // point number within the patch system. @@ -46,55 +69,49 @@ // very general, but matches our present use in this thorn. // +//****************************************************************************** + +// ABC to define Jacobian-traversal interface class Jacobian { public: - // access a matrix element via row/column indices - virtual const fp& operator()(int II, int JJ) const = 0; // rvalue - virtual fp& operator()(int II, int JJ) = 0; // lvalue + // basic meta-info + patch_system& my_patch_system() const { return ps_; } + int NN() const { return NN_; } - // access a matrix element via row index and column (patch,irho,isigma) - const fp& operator() // rvalue - (int II, const patch& JJ_patch, int JJ_irho, int JJ_isigma) + // convert (patch,irho,isigma) <--> row/column index + int II_of_patch_irho_isigma(const patch& p, int irho, int isigma) const - { - const int JJ - = ps_.gpn_of_patch_irho_isigma(JJ_patch, JJ_irho,JJ_isigma); - return operator()(II,JJ); - } - fp& operator() // lvalue - (int II, const patch& JJ_patch, int JJ_irho, int JJ_isigma) - { - const int JJ - = ps_.gpn_of_patch_irho_isigma(JJ_patch, JJ_irho,JJ_isigma); - return operator()(II,JJ); - } - - // is a given element explicitly stored, or implicitly 0 via sparsity - virtual bool is_explicitly_stored(int II, int JJ) const = 0; + { return ps_.gpn_of_patch_irho_isigma(p, irho,isigma); } + const patch& patch_irho_isigma_of_II(int II, int& irho, int& isigma) + const + { return ps_.patch_irho_isigma_of_gpn(II, irho,isigma); } - // zero the whole matrix or a single row + // zero the entire matrix virtual void zero_matrix() = 0; - virtual void zero_row(int II) = 0; - // solve linear system J.x = rhs - // ... rhs and x are nominal-grid gridfns - // ... may modify Jacobian matrix (eg for LU decomposition) - // ... returns 0.0 if matrix is numerically singular - // condition number (> 0.0) if known - // -1.0 if condition number is unknown - virtual fp solve_linear_system(int rhs_gfn, int x_gfn) = 0; + // set a matrix element to a specified value + virtual void set_element(int II, int JJ, fp value) = 0; - // basic meta-info - patch_system& my_patch_system() const { return ps_; } - int NN() const { return NN_; } + // sum a value into a matrix element + virtual void sum_into_element(int II, int JJ, fp value) = 0; +protected: // constructor, destructor - // ... for analyze-factor-solve sparse-matrix derived classes, - // construction may include the "analyze" operation - Jacobian(patch_system& ps); + Jacobian(patch_system& ps) + : ps_(ps), + NN_(ps.N_grid_points()) + { } +public: virtual ~Jacobian() { } +private: + // we forbid copying and passing by value + // by declaring the copy constructor and assignment operator + // private, but never defining them + Jacobian(const Jacobian &rhs); + Jacobian& operator=(const Jacobian &rhs); + protected: patch_system& ps_; int NN_; // ps_.N_grid_points() @@ -102,65 +119,80 @@ protected: //****************************************************************************** -#ifdef HAVE_DENSE_JACOBIAN -// -// This class stores the Jacobian as a dense matrix in Fortran (column) -// order. -// -class dense_Jacobian +// ABC to accumulate Jacobian sparsity info +class Jacobian_sparsity : public Jacobian { public: - // access a matrix element via row/column indices - const fp& operator()(int II, int JJ) const // rvalue - { return matrix_(JJ,II); } - fp& operator()(int II, int JJ) // lvalue - { return matrix_(JJ,II); } - - bool is_explicitly_stored(int II, int JJ) const - { return true; } - - // zero the whole matrix or a single row - void zero_matrix(); - void zero_row(int II); + Jacobian_sparsity(patch_system& ps, + bool print_msg_flag = false) + : Jacobian(ps) + { } + virtual ~Jacobian_sparsity() { } + }; - // solve linear system J.x = rhs via LAPACK - // ... rhs and x are nominal-grid gridfns - // ... overwrites Jacobian matrix with its LU decomposition - // ... returns 0.0 if matrix is numerically singular, or - // condition number (> 0.0) - fp solve_linear_system(int rhs_gfn, int x_gfn); +//****************************************************************************** - // constructor, destructor +// Jacobian_matrix -- ABC to store/use Jacobian matrix +class Jacobian_matrix + : public Jacobian + { public: - dense_Jacobian(patch_system& ps); - ~dense_Jacobian(); + // access (get) a matrix element + virtual fp element(int II, int JJ) const = 0; -private: - // subscripts are (JJ,II) (both 0-origin) - jtutil::array2d<fp> matrix_; + // is a given element explicitly stored, or implicitly 0 via sparsity + virtual bool is_explicitly_stored(int II, int JJ) const = 0; - // pivot vector for LAPACK routines - integer *pivot_; + // solve linear system J.x = rhs + // ... rhs and x are nominal-grid gridfns + // ... may modify Jacobian matrix (eg for LU decomposition) + // ... returns 0.0 if matrix is numerically singular + // condition number (> 0.0) if known + // -1.0 if condition number is unknown + virtual fp solve_linear_system(int rhs_gfn, int x_gfn) = 0; - // work vectors for LAPACK condition number computation - integer *iwork_; - fp *rwork_; + // constructor, destructor + // ... note Jacobian_sparsity argument of constructor is non-const + // to allow us to take over ownership of data structures + Jacobian_matrix(patch_system& ps, + Jacobian_sparsity& JS, + bool print_msg_flag = false) + : Jacobian(ps) + { } + virtual ~Jacobian_matrix() { } }; -#endif // HAVE_DENSE_JACOBIAN //****************************************************************************** +//****************************************************************************** +//****************************************************************************** enum Jacobian_type { #ifdef HAVE_DENSE_JACOBIAN - Jacobian_type__dense_matrix // no comma on last entry in enum + Jacobian_type__dense_matrix, +#endif +#ifdef HAVE_ROW_SPARSE_JACOBIAN + Jacobian_type__row_sparse_matrix // no comma on last entry in enum #endif }; +// +// prototypes of various Jacobian-related functions +// + // decode string into internal enum enum Jacobian_type decode_Jacobian_type(const char Jacobian_type_string[]); -// construct and return new-allocated Jacobian object of specified type -Jacobian& new_Jacobian(patch_system& ps, enum Jacobian_type Jac_type); +// "object factory" routines to construct and return +// pointers to new-allocated objects of specified derived types +class Jacobian_sparsity; +class Jacobian_matrix; +Jacobian_sparsity* new_Jacobian_sparsity(enum Jacobian_type Jac_type, + patch_system& ps, + bool print_msg_flag = false); +Jacobian_matrix* new_Jacobian_matrix(enum Jacobian_type Jac_type, + patch_system& ps, + Jacobian_sparsity& JS, + bool print_msg_flag = false); diff --git a/src/elliptic/make.code.defn b/src/elliptic/make.code.defn index 3bc82f4..28308a0 100644 --- a/src/elliptic/make.code.defn +++ b/src/elliptic/make.code.defn @@ -2,7 +2,9 @@ # $Header$ # Source files in this directory -SRCS = Jacobian.cc gecon_wrapper.F77 +SRCS = Jacobian.cc \ + dense_Jacobian.cc gecon_wrapper.F77 \ + row_sparse_Jacobian.cc # Subdirectories containing source files SUBDIRS = diff --git a/src/gr/expansion_Jacobian.cc b/src/gr/expansion_Jacobian.cc index a843e4e..5faf059 100644 --- a/src/gr/expansion_Jacobian.cc +++ b/src/gr/expansion_Jacobian.cc @@ -94,7 +94,7 @@ bool expansion_Jacobian(patch_system& ps, const struct cactus_grid_info& cgi, const struct geometry_info& gi, const struct Jacobian_info& Jacobian_info, - bool print_msg_flag) + bool print_msg_flag /* = false */) { switch (Jacobian_info.Jacobian_method) { @@ -136,8 +136,8 @@ return true; // *** NORMAL RETURN *** // // This function computes the Jacobian matrix of the expansion Theta(h) -// by numerical perturbation of the Theta(h) function. The algorithm is as -// follows: +// by numerical perturbation of the Theta(h) function. The algorithm is +// as follows: // // we assume that Theta = Theta(h) has already been evaluated // save_Theta = Theta @@ -161,6 +161,7 @@ return true; // *** NORMAL RETURN *** // // Outputs: // The Jacobian matrix is stored in the Jacobian object Jac. +// As implied by the above algorithm, it's traversed by columns. // // Results: // This function returns true for a successful computation, or false @@ -219,7 +220,7 @@ ps.gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta); x_irho,x_isigma); const fp new_Theta = xp.gridfn(gfns::gfn__Theta, x_irho,x_isigma); - Jac(II,JJ) = (new_Theta - old_Theta) / epsilon; + Jac.set_element(II,JJ, (new_Theta - old_Theta) / epsilon); } } } @@ -238,9 +239,8 @@ return true; // *** NORMAL RETURN *** // // This function computes the partial derivative terms in the Jacobian // matrix of the expansion Theta(h), by symbolic differentiation from -// the Jacobian coefficient (angular) gridfns. The computation is done a -// Jacobian row at a time, using equation (25) of my 1996 apparent horizon -// finding paper. +// the Jacobian coefficient (angular) gridfns. The Jacobian is traversed +// by rows, using equation (25) of my 1996 apparent horizon finding paper. // // Inputs (angular gridfns, on ghosted grid): // h # shape of trial surface @@ -316,7 +316,11 @@ ps.compute_synchronize_Jacobian(); * xp.partial_rho_rho_coeff(m_irho); const fp Jac_sum = Jac_rho + Jac_rho_rho; if (xp.is_in_nominal_grid(xm_irho, x_isigma)) - then Jac(II, xp,xm_irho,x_isigma) += Jac_sum; + then { + const int xm_JJ + = Jac.II_of_patch_irho_isigma(xp,xm_irho,x_isigma); + Jac.sum_into_element(II, xm_JJ, Jac_sum); + } else add_ghost_zone_Jacobian (ps, Jac, Jac_sum, @@ -338,7 +342,11 @@ ps.compute_synchronize_Jacobian(); * xp.partial_sigma_sigma_coeff(m_isigma); const fp Jac_sum = Jac_sigma + Jac_sigma_sigma; if (xp.is_in_nominal_grid(x_irho, xm_isigma)) - then Jac(II, xp,x_irho,xm_isigma) += Jac_sum; + then { + const int xm_JJ + = Jac.II_of_patch_irho_isigma(xp, x_irho, xm_isigma); + Jac.sum_into_element(II, xm_JJ, Jac_sum); + } else add_ghost_zone_Jacobian (ps, Jac, Jac_sum, @@ -363,7 +371,11 @@ ps.compute_synchronize_Jacobian(); = Jacobian_coeff_rho_sigma * xp.partial_rho_sigma_coeff(m_irho, m_isigma); if (xp.is_in_nominal_grid(xm_irho, xm_isigma)) - then Jac(II, xp,xm_irho,xm_isigma) += Jac_rho_sigma; + then { + const int xm_JJ + = Jac.II_of_patch_irho_isigma(xp, xm_irho, xm_isigma); + Jac.sum_into_element(II, xm_JJ, Jac_rho_sigma); + } else { const ghost_zone& xmgz = xp.corner_ghost_zone_containing_point @@ -428,7 +440,7 @@ const patch_edge& ye = ps.synchronize_Jacobian(xmgz, y_iperp, y_posn, min_ym, max_ym, Jacobian_buffer); -const patch& yp = ye.my_patch(); +patch& yp = ye.my_patch(); // add the Jacobian contributions from the ym points for (int ym = min_ym ; ym <= max_ym ; ++ym) @@ -436,8 +448,8 @@ const patch& yp = ye.my_patch(); const int y_ipar = y_posn + ym; const int y_irho = ye. irho_of_iperp_ipar(y_iperp,y_ipar); const int y_isigma = ye.isigma_of_iperp_ipar(y_iperp,y_ipar); - const int JJ = ps.gpn_of_patch_irho_isigma(yp, y_irho, y_isigma); - Jac(x_II,JJ) += mol*Jacobian_buffer(ym); + const int y_JJ = Jac.II_of_patch_irho_isigma(yp, y_irho, y_isigma); + Jac.sum_into_element(x_II, y_JJ, mol*Jacobian_buffer(ym)); } } } @@ -498,7 +510,7 @@ if (! expansion(ps, cgi, gi)) irho,isigma); const fp new_Theta = p.gridfn(gfns::gfn__Theta, irho,isigma); - Jac(II,II) += (new_Theta - old_Theta) / epsilon; + Jac.sum_into_element(II,II, (new_Theta - old_Theta) / epsilon); } } } diff --git a/src/gr/gr.hh b/src/gr/gr.hh index b5e06f7..43df6ae 100644 --- a/src/gr/gr.hh +++ b/src/gr/gr.hh @@ -170,7 +170,7 @@ bool expansion_Jacobian(patch_system& ps, const struct cactus_grid_info& cgi, const struct geometry_info& gi, const struct Jacobian_info& Jacobian_info, - bool print_msg_flag); + bool print_msg_flag = false); // Schwarzschild_EF.cc void Schwarzschild_EF_geometry(patch_system& ps, @@ -183,3 +183,9 @@ enum geometry_method bool geometry_method_is_interp(enum geometry_method geometry_method); enum Jacobian_method decode_Jacobian_method(const char Jacobian_method_string[]); +Jacobian_matrix* + create_Jacobian_matrix(patch_system& ps, + const struct cactus_grid_info& cgi, + const struct geometry_info& gi, + const struct Jacobian_info& Jac_info, + bool print_msg_flag = false); diff --git a/src/gr/misc.cc b/src/gr/misc.cc index 28cdfbd..3951331 100644 --- a/src/gr/misc.cc +++ b/src/gr/misc.cc @@ -4,6 +4,7 @@ // decode_geometry_method - decode the geometry_method parameter // geomery_method_is_interp - does enum geometry_method specify interpolation? // decode_geometry_method - decode the geometry_method parameter +// create_Jacobian_matrix - create a new Jacobian_matrix object // #include <stdio.h> @@ -103,3 +104,49 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, "decode_Jacobian_method(): unknown Jacobian_method_string=\"%s\"!", Jacobian_method_string); /*NOTREACHED*/ } + +//****************************************************************************** + +// +// This function creates a new-allocated Jacobian matrix and returns a +// pointer to it. In more detail, it constructs a Jacobian_sparsity object, +// calls expansion() and then expansion_Jacobian() to traverse the Jacobian +// matrix (recording the Jacobian's sparsity pattern in the process), then +// uses this sparsity pattern to construct the Jacobian matrix itself +// (and destroys the Jacobian_sparsity object). +// +Jacobian_matrix* + create_Jacobian_matrix(patch_system& ps, + const struct cactus_grid_info& cgi, + const struct geometry_info& gi, + const struct Jacobian_info& Jac_info, + bool print_msg_flag /* = false */) +{ +if (print_msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, + "traversing Jacobian matrix (recording sparsity pattern)"); +Jacobian_sparsity* JS_ptr + = new_Jacobian_sparsity(Jac_info.Jacobian_storage_method, + ps, + print_msg_flag); + +// it's safest to compute expansion() before doing the traversal +// to make sure all the Jacobian coefficients are set to sane values +expansion(ps, + cgi, gi, + true); // yes, compute Jacobian coeffs + +// now do the actual traversal (==> records the sparsity pattern) +expansion_Jacobian(ps, *JS_ptr, + cgi, gi, Jac_info); + +if (print_msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, + "constructing Jacobian matrix"); +Jacobian_matrix* Jac_ptr + = new_Jacobian_matrix(Jac_info.Jacobian_storage_method, + ps, *JS_ptr, + print_msg_flag); +delete JS_ptr; +return Jac_ptr; +} diff --git a/src/include/config.h b/src/include/config.h index 511081f..492da5a 100644 --- a/src/include/config.h +++ b/src/include/config.h @@ -61,7 +61,9 @@ typedef CCTK_INT integer; /* * What types of Jacobian matrix storage do we want to compile in * support for? N.b. each of these requires linking with the corresponding - * linear-solver library; see "../make.configuration.defn" for details - * on these libraries. + * linear-solver library; see "../elliptic/Jacobian.hh" for details + * on the storage formats, or "../make.configuration.defn" (if it exists) + * for details on the libraries. */ #define HAVE_DENSE_JACOBIAN +#define HAVE_ROW_SPARSE_JACOBIAN |