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, 27 insertions, 22 deletions

diff --git a/src/gr/Newton.cc b/src/gr/Newton.cc
index 011f34c..4c2eb04 100644
--- a/src/gr/Newton.cc
+++ b/src/gr/Newton.cc
@@ -51,44 +51,39 @@ namespace {
 
 //
 // This function solves H(h) = 0 for h via Newton's method.
-// At the moment the NP Jacobian is used.
 //
-void Newton_solve(patch_system& ps,
+bool Newton_solve(patch_system& ps,
 		  const struct cactus_grid_info& cgi,
 		  const struct geometry_interpolator_info& gii,
-		  const char Jacobian_type[],
-		  fp perturbation_amplitude,
-		  int max_Newton_iterations,
-		  fp H_norm_for_convergence)
+		  const struct Jacobian_info& Jac_info,
+		  const struct solver_info& solver_info,
+		  Jacobian& Jac)
 {
-Jacobian& Jac = create_Jacobian(ps, Jacobian_type);
-
 	for (int iteration = 1 ;
-	     iteration <= max_Newton_iterations ;
+	     iteration <= solver_info.max_Newton_iterations ;
 	     ++iteration)
 	{
 	CCTK_VInfo(CCTK_THORNSTRING,
 		   "Newton iteration %d", iteration);
 
 	jtutil::norm<fp> H_norms;
-	horizon_function(ps, cgi, gii, true, H_norms);
+	horizon_function(ps, cgi, gii, true, &H_norms);
 	CCTK_VInfo(CCTK_THORNSTRING,
 		   "   H rms-norm=%.2e, infinity-norm=%.2e",
 		   H_norms.rms_norm(), H_norms.infinity_norm());
 
-	if (H_norms.infinity_norm() <= H_norm_for_convergence)
+	if (H_norms.infinity_norm() <= solver_info.H_norm_for_convergence)
 	   then {
-		CCTK_VInfo(CCTK_THORNSTRING, "==> finished!");
-		return;					// *** NORMAL RETURN ***
+		CCTK_VInfo(CCTK_THORNSTRING,
+			   "==> finished (||H|| < %g)",
+			   double(solver_info.H_norm_for_convergence));
+		return true;				// *** NORMAL RETURN ***
 		}
 
 	// take a Newton step
-	horizon_Jacobian_NP(ps, cgi, gii, Jac, perturbation_amplitude);
-	jtutil::array_scale(Jac.NN(),
-			    -1.0,
-			    ps.gridfn_data(nominal_gfns::gfn__H));
+	horizon_Jacobian(ps, cgi, gii, Jac_info, Jac);
 	CCTK_VInfo(CCTK_THORNSTRING,
-		   "solving linear system");
+		   "solving linear system for %d unknowns", Jac.NN());
 	Jac.solve_linear_system(nominal_gfns::gfn__H,		// rhs gridfn
 				nominal_gfns::gfn__Delta_h);	// soln gridfn
 	CCTK_VInfo(CCTK_THORNSTRING,
@@ -108,16 +103,26 @@ Jacobian& Jac = create_Jacobian(ps, Jacobian_type);
 					    irho,isigma);
 		Delta_h_norms.data(Delta_h);
 
-		p.ghosted_gridfn(ghosted_gfns::gfn__h, irho,isigma) += Delta_h;
+		p.ghosted_gridfn(ghosted_gfns::gfn__h, irho,isigma) -= Delta_h;
 		}
 		}
 		}
 	CCTK_VInfo(CCTK_THORNSTRING,
 		   "moved h by Delta_h (rhs-norm=%.2e, infinity-norm=%.2e)",
 		   Delta_h_norms.rms_norm(), Delta_h_norms.infinity_norm());
+
+	if (Delta_h_norms.infinity_norm() <= solver_info
+					     .Delta_h_norm_for_convergence)
+	   then {
+		CCTK_VInfo(CCTK_THORNSTRING,
+			   "==> finished (||Delta_h|| < %g)",
+			   double(solver_info.Delta_h_norm_for_convergence));
+		return true;				// *** NORMAL RETURN ***
+		}
 	}
 
-CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
-	   "Newton_solve: failed to converge in %d iterations!",
-	   max_Newton_iterations);
+CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+	   "Newton_solve: no convergence in %d iterations!\n",
+	   solver_info.max_Newton_iterations);
+return false;						// *** ERROR RETURN ***
 }