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

1 files changed, 11 insertions, 135 deletions

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);
-	  };
-
 //******************************************************************************
 //******************************************************************************
 //******************************************************************************
@@ -100,130 +91,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.
 //
 dense_Jacobian::dense_Jacobian(patch_system& ps)
@@ -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;