From b3c9f3f1269b213537b25a1ff0123871a5e07f9f Mon Sep 17 00:00:00 2001 From: jthorn Date: Mon, 22 Jul 2002 12:29:42 +0000 Subject: add a whole bunch of changes from working at home --> AHFinderDirect now finds AH correctly for Kerr/offset!!! git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@648 f88db872-0e4f-0410-b76b-b9085cfa78c5 --- src/elliptic/Jacobian.cc | 238 +++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 231 insertions(+), 7 deletions(-) (limited to 'src/elliptic/Jacobian.cc') diff --git a/src/elliptic/Jacobian.cc b/src/elliptic/Jacobian.cc index 9b3e4d8..4d3e411 100644 --- a/src/elliptic/Jacobian.cc +++ b/src/elliptic/Jacobian.cc @@ -4,12 +4,22 @@ // // Jacobian::Jacobian // +// print_Jacobian +// print_Jacobians +// // dense_Jacobian::dense_Jacobian +// dense_Jacobian::~dense_Jacobian // dense_Jacobian::zero // dense_Jacobian::zero_row +// dense_Jacobian::solve_linear_system // #include +using std::fopen; +using std::printf; +using std::fprintf; +using std::fclose; +using std::FILE; #include #include #include @@ -36,6 +46,42 @@ using jtutil::error_exit; #include "../util/patch_system.hh" #include "Jacobian.hh" +//#include "lapack.h" +//************************************** +/* lapack.h -- C/C++ prototypes for (some) LAPACK routines */ +/* $Id$ */ + +/* + * prerequisites: + * "cctk.h" + * "config.hh" // for "integer" = Fortran integer + */ + +#ifdef __cplusplus +extern "C" { +#endif + +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); + +#ifdef __cplusplus + }; /* extern "C" */ +#endif +//************************************** + +namespace { +void print_Jacobians_internal(const char file_name[], + const Jacobian& SD_Jac, const Jacobian& NP_Jac, + bool pair_flag); + }; //****************************************************************************** //****************************************************************************** @@ -44,31 +90,166 @@ using jtutil::error_exit; // // This function constructs a Jacobian object. // -Jacobian::Jacobian(const patch_system& ps) - : ps_(ps) +Jacobian::Jacobian(patch_system& ps) + : ps_(ps), + NN_(ps.N_grid_points()) { } //****************************************************************************** //****************************************************************************** //****************************************************************************** +// +// This function prints a Jacobian matrix to a named output file. +// +void print_Jacobian(const char file_name[], const Jacobian& Jac) +{ +print_Jacobians_internal(file_name, Jac, Jac, false); +} + +//****************************************************************************** + +// +// This function prints a pair of Jacobian matrices (and their difference) +// to a named output file. +// +void print_Jacobians(const char file_name[], + const Jacobian& SD_Jac, const Jacobian& NP_Jac) +{ +print_Jacobians_internal(file_name, SD_Jac, NP_Jac, true); +} + +//****************************************************************************** + +// +// If pair_flag = false, this prints SD_Jac. +// If pair_flag = true, this prints both Jacobians, and the error in SD_Jac. +// +namespace { +void print_Jacobians_internal(const char file_name[], + const Jacobian& SD_Jac, const Jacobian& NP_Jac, + bool pair_flag) +{ +const patch_system& ps = SD_Jac.my_patch_system(); + +FILE *fileptr = fopen(file_name, "w"); +if (fileptr == NULL) + then error_exit(ERROR_EXIT, +"***** dense_Jacobian::print(): can't open output file \"%s\"!", + file_name); /*NOTREACHED*/ + +fprintf(fileptr, "# column 1 = x II\n"); +fprintf(fileptr, "# column 2 = x patch number\n"); +fprintf(fileptr, "# column 3 = x irho\n"); +fprintf(fileptr, "# column 4 = x isigma\n"); +fprintf(fileptr, "# column 5 = y JJ\n"); +fprintf(fileptr, "# column 6 = y patch number\n"); +fprintf(fileptr, "# column 7 = y irho\n"); +fprintf(fileptr, "# column 8 = y isigma\n"); +if (pair_flag) + then { + fprintf(fileptr, "# column 9 = SD_Jac(II,JJ)\n"); + fprintf(fileptr, "# column 10 = NP_Jac(II,JJ)\n"); + fprintf(fileptr, "# column 11 = abs error\n"); + fprintf(fileptr, "# column 12 = rel error\n"); + } + else fprintf(fileptr, "# column 9 = Jac(II,JJ)\n"); + + for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn) + { + patch& xp = ps.ith_patch(xpn); + + for (int x_irho = xp.min_irho() ; x_irho <= xp.max_irho() ; ++x_irho) + { + for (int x_isigma = xp.min_isigma() ; + x_isigma <= xp.max_isigma() ; + ++x_isigma) + { + const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho,x_isigma); + + for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn) + { + patch& yp = ps.ith_patch(ypn); + + for (int y_irho = yp.min_irho() ; + y_irho <= yp.max_irho() ; + ++y_irho) + { + for (int y_isigma = yp.min_isigma() ; + y_isigma <= yp.max_isigma() ; + ++y_isigma) + { + const int JJ = ps.gpn_of_patch_irho_isigma(yp, y_irho,y_isigma); + + if (! SD_Jac.is_explicitly_stored(II,JJ)) + then continue; // skip sparse points + + const fp SD = SD_Jac(II,JJ); + const fp NP = NP_Jac(II,JJ); + const fp abs_error = SD - NP; + + if (pair_flag ? ((SD == 0.0) && (NP == 0.0)) + : (SD == 0.0)) + then continue; // skip zero values (if == ) + + const fp abs_SD = jtutil::abs(SD); + const fp abs_NP = jtutil::abs(NP); + const fp rel_error = abs_error / jtutil::max(abs_SD, abs_NP); + + fprintf(fileptr, + "%d %d %d %d\t%d %d %d %d\t", + II, xpn, x_irho, x_isigma, + JJ, ypn, y_irho, y_isigma); + if (pair_flag) + then fprintf(fileptr, + "%.10g\t%.10g\t%e\t%e\n", + double(SD), double(NP), + double(abs_error), double(rel_error)); + else fprintf(fileptr, + "%.10g\n", + double(SD)); + } + } + } + } + } + } + +fclose(fileptr); +} + }; + +//****************************************************************************** +//****************************************************************************** +//****************************************************************************** + // // This function constructs a dense_Jacobian object. // -dense_Jacobian::dense_Jacobian(const patch_system& ps) +dense_Jacobian::dense_Jacobian(patch_system& ps) : Jacobian(ps), - matrix_(0, ps.N_grid_points(), - 0, ps.N_grid_points()) + matrix_(0,NN()-1, 0,NN()-1), + pivot_(new integer[NN()]) { } //****************************************************************************** +// +// THis function destroys a dense_Jacobian object. +// +dense_Jacobian::~dense_Jacobian() +{ +delete[] pivot_; +} + +//****************************************************************************** + // // This function zeros a dense_Jacobian object. // void dense_Jacobian::zero() { -// FIXME FIXME +jtutil::array_zero(matrix_.N_array(), matrix_.data_array()); } //****************************************************************************** @@ -78,5 +259,48 @@ void dense_Jacobian::zero() // void dense_Jacobian::zero_row(int II) { -// FIXME FIXME + for (int JJ = 0 ; JJ < NN() ; ++JJ) + { + matrix_(JJ,II) = 0.0; + } +} + +//****************************************************************************** + +// +// This function solves the linear system J.x = rhs, with rhs and x +// being nominal-grid gridfns. The computation is done using the LAPACK +// [sd]gesv() subroutines; which modify the Jacobian matrix J in-place +// for the LU decomposition. +// +void 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); + +// +// [sd]gesv() use an "in out" design, where the same argument is used for +// both rhs and x ==> first copy rhs to x so we can pass that to [sd]gesv() +// +jtutil::array_copy(NN(), rhs, x); + +integer N = NN(); +integer NRHS = 1; +integer info; + +#ifdef FP_IS_FLOAT + CCTK_FNAME(sgesv)(&N, &NRHS, matrix_.data_array(), &N, pivot_, x, &N, &info); +#endif +#ifdef FP_IS_DOUBLE + CCTK_FNAME(dgesv)(&N, &NRHS, matrix_.data_array(), &N, pivot_, x, &N, &info); +#endif + +if (info != 0) + then error_exit(ERROR_EXIT, +"\n" +"***** dense_Jacobian::solve_linear_system(rhs_gfn=%d, x_gfn=%d):\n" +" error return info=%d from [sd]gesv() LAPACK routine!\n" +, + rhs_gfn, x_gfn, + int(info)); /*NOTREACHED*/ } -- cgit v1.2.3