// Newton.cc -- solve H(h) = 0 via Newton's method
// $Id$
//
// <<<prototypes for functions local to this file>>>
// Newton_solve - driver to solve H(h) = 0 via Newton's method
//

#include <stdio.h>
#include <assert.h>
#include <math.h>
#include <vector>

#include "util_Table.h"
#include "cctk.h"
#include "cctk_Arguments.h"

#include "stdc.h"
#include "config.hh"
#include "../jtutil/util.hh"
#include "../jtutil/array.hh"
#include "../jtutil/cpm_map.hh"
#include "../jtutil/linear_map.hh"
using jtutil::error_exit;

#include "../util/coords.hh"
#include "../util/grid.hh"
#include "../util/fd_grid.hh"
#include "../util/patch.hh"
#include "../util/patch_edge.hh"
#include "../util/patch_interp.hh"
#include "../util/ghost_zone.hh"
#include "../util/patch_system.hh"

#include "../elliptic/Jacobian.hh"

#include "gfns.hh"
#include "AHFinderDirect.hh"

//******************************************************************************

//
// ***** prototypes for functions local to this file *****
//

namespace {
	  }

//******************************************************************************
//******************************************************************************
//******************************************************************************

//
// This function solves H(h) = 0 for h via Newton's method.
//
bool Newton_solve(patch_system& ps,
		  const struct cactus_grid_info& cgi,
		  const struct geometry_info& gi,
		  const struct Jacobian_info& Jacobian_info,
		  const struct solver_info& solver_info,
		  Jacobian& Jac)
{
	for (int iteration = 1 ;
	     iteration <= solver_info.max_Newton_iterations ;
	     ++iteration)
	{
	CCTK_VInfo(CCTK_THORNSTRING,
		   "Newton iteration %d", iteration);

	jtutil::norm<fp> H_norms;
	horizon_function(ps, cgi, gi, true, true, &H_norms);
	CCTK_VInfo(CCTK_THORNSTRING,
		   "   H rms-norm=%.2e, infinity-norm=%.2e",
		   H_norms.rms_norm(), H_norms.infinity_norm());

	if (solver_info.output_h_and_H_at_each_Newton_iteration)
	   then {
		const char N_buffer = 100;
		char file_name_buffer[100];

		snprintf(file_name_buffer, N_buffer,
			 "%s.it%d",
			 solver_info.h_file_name, iteration);
		CCTK_VInfo(CCTK_THORNSTRING,
			   "   writing h to \"%s\"",
			   file_name_buffer);
		ps.print_ghosted_gridfn_with_xyz(gfns::gfn__h,
						 true, gfns::gfn__h,
						 file_name_buffer,
						 false);   // no ghost zones
        
		snprintf(file_name_buffer, N_buffer,
			 "%s.it%d",
			 solver_info.H_of_h_file_name, iteration);
		CCTK_VInfo(CCTK_THORNSTRING,
			   "   writing H(h) to \"%s\"",
			   file_name_buffer);
		ps.print_gridfn_with_xyz(gfns::gfn__H,
					 true, gfns::gfn__h,
					 file_name_buffer);
		}

	if (H_norms.infinity_norm() <= solver_info.H_norm_for_convergence)
	   then {
		CCTK_VInfo(CCTK_THORNSTRING,
			   "==> finished (||H|| < %g)",
			   double(solver_info.H_norm_for_convergence));
		return true;				// *** NORMAL RETURN ***
		}

	// compute the Newton step
	horizon_Jacobian(ps, cgi, gi, Jacobian_info, Jac);
	CCTK_VInfo(CCTK_THORNSTRING,
		   "solving linear system for Delta_h (%d unknowns)",
		   Jac.NN());
	const fp rcond = Jac.solve_linear_system(gfns::gfn__H,
						 gfns::gfn__Delta_h);
	if	(rcond == 0.0)
	   then {
		CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
		   "Newton_solve: Jacobian matrix is numerically singular!");
		return false;				// *** ERROR RETURN ***
		}
	if (rcond > 0.0)
	   then CCTK_VInfo(CCTK_THORNSTRING,
			   "done solving linear system (condition number %.2e)",
			   double(rcond));
	   else CCTK_VInfo(CCTK_THORNSTRING,
			   "done solving linear system");

	// take the Newton step
	jtutil::norm<fp> Delta_h_norms;
		for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
		{
		patch& p = ps.ith_patch(pn);

		for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho)
		{
		for (int isigma = p.min_isigma() ;
		     isigma <= p.max_isigma() ;
		     ++isigma)
		{
		const fp Delta_h = p.gridfn(gfns::gfn__Delta_h, irho,isigma);
		Delta_h_norms.data(Delta_h);

		p.ghosted_gridfn(gfns::gfn__h, irho,isigma) -= Delta_h;
		}
		}
		}
	CCTK_VInfo(CCTK_THORNSTRING,
		   "h += Delta_h (rms-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));
		if (solver_info.final_H_update_if_exit_x_H_small)
		   then {
			CCTK_VInfo(CCTK_THORNSTRING,
				   "    doing final H(h) evaluation");
			horizon_function(ps, cgi, gi, false, true);
			}
		return true;				// *** NORMAL RETURN ***
		}
	}


CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
	   "Newton_solve: no convergence in %d iterations!",
	   solver_info.max_Newton_iterations);
if (solver_info.final_H_update_if_exit_x_H_small)
   then {
	CCTK_VInfo(CCTK_THORNSTRING,
		   "    doing final H(h) evaluation");
	horizon_function(ps, cgi, gi, false, true);
	}
return false;						// *** ERROR RETURN ***
}