From 73d9cdfe58e99305ede9c394f72712e7f759b5a2 Mon Sep 17 00:00:00 2001 From: jthorn Date: Mon, 22 Jul 2002 17:33:21 +0000 Subject: various changes including d/dr terms in Jacobian by numerical perturbation, tweak I/O parameters, move printing Jacobian out of Jacobian class into test driver, drop unused array BLAS routines in jtutil:: git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@654 f88db872-0e4f-0410-b76b-b9085cfa78c5 --- src/elliptic/Jacobian.cc | 146 ++++------------------------------------------- 1 file changed, 11 insertions(+), 135 deletions(-) (limited to 'src/elliptic/Jacobian.cc') diff --git a/src/elliptic/Jacobian.cc b/src/elliptic/Jacobian.cc index 4d3e411..58009ec 100644 --- a/src/elliptic/Jacobian.cc +++ b/src/elliptic/Jacobian.cc @@ -4,9 +4,6 @@ // // Jacobian::Jacobian // -// print_Jacobian -// print_Jacobians -// // dense_Jacobian::dense_Jacobian // dense_Jacobian::~dense_Jacobian // dense_Jacobian::zero @@ -77,12 +74,6 @@ void CCTK_FCALL #endif //************************************** -namespace { -void print_Jacobians_internal(const char file_name[], - const Jacobian& SD_Jac, const Jacobian& NP_Jac, - bool pair_flag); - }; - //****************************************************************************** //****************************************************************************** //****************************************************************************** @@ -99,130 +90,6 @@ Jacobian::Jacobian(patch_system& ps) //****************************************************************************** //****************************************************************************** -// -// 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. // @@ -249,7 +116,13 @@ delete[] pivot_; // void dense_Jacobian::zero() { -jtutil::array_zero(matrix_.N_array(), matrix_.data_array()); + for (int JJ = 0 ; JJ < NN() ; ++JJ) + { + for (int II = 0 ; II < NN() ; ++II) + { + matrix_(JJ,II) = 0.0; + } + } } //****************************************************************************** @@ -282,7 +155,10 @@ 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); + for (int II = 0 ; II < NN() ; ++II) + { + x[II] = rhs[II]; + } integer N = NN(); integer NRHS = 1; -- cgit v1.2.3