//#####################################################################
// Copyright 2004-2008, Ron Fedkiw, Geoffrey Irving, Michael Lentine, Frank Losasso, Andrew Selle.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class PCG_SPARSE
//#####################################################################
//
// Solve Ax=b with Preconditioned Conjugate Gradient, Modified Incomplete Cholesky preconditioner.
// In the Neumann case, the input A, b, & x need to sum to 0 for compatibility.
// Input A as a sparse (symmetric positive definite) matrix & b as a vector. Input x as a vector. 
// Input tolerance for stopping iteration when the maximum residual is less than or equal to the tolerance, i.e. max|r|=max|b-Ax| <= tolerance. Otherwise stops at max iterations.
//               
//#####################################################################
#ifndef __PCG_SPARSE__
#define __PCG_SPARSE__

#include <Core/Krylov_Solvers/KRYLOV_SOLVER.h>
namespace PhysBAM{

template<class T> class SPARSE_MATRIX_FLAT_MXN;

template<class T>
class PCG_SPARSE
{
public:
    bool enforce_compatibility; // for Neumann boundary conditions
    bool remove_null_space_solution_component; 
    bool show_residual,show_results,ensure_convergence;
    bool incomplete_cholesky; // true when using this preconditioner or the modified version below
    bool modified_incomplete_cholesky; // true when using this preconditioner
    T modified_incomplete_cholesky_coefficient;
    int maximum_iterations;
    int cg_restart_iterations;
    T preconditioner_zero_tolerance,preconditioner_zero_replacement;
    KRYLOV_SOLVER_TYPE evolution_solver_type;

    PCG_SPARSE()
        :cg_restart_iterations(0),preconditioner_zero_tolerance((T)1e-8),preconditioner_zero_replacement((T)1e-8),evolution_solver_type(krylov_solver_cg)
    {
        Enforce_Compatibility(false);
        Remove_Null_Space_Solution_Component(false);
        Show_Residuals(false);Show_Results(false);
        ensure_convergence=false;
        Use_Modified_Incomplete_Cholesky();
        Set_Maximum_Iterations();
    }

    PCG_SPARSE(const PCG_SPARSE&) = delete;
    void operator=(const PCG_SPARSE&) = delete;
    
    ~PCG_SPARSE()
    {}

    void Enforce_Compatibility(const bool enforce_compatibility_input=true) //  A, b, & x inputs need to sum to 0 for compatibility
    {enforce_compatibility=enforce_compatibility_input;}

    void Remove_Null_Space_Solution_Component(const bool remove_null_space_solution_component_input=true)
    {remove_null_space_solution_component=remove_null_space_solution_component_input;}
    
    void Show_Residuals(const bool show_residual_input=true)
    {show_residual=show_residual_input;}

    void Show_Results(const bool show_results_input=true)
    {show_results=show_results_input;}

    void Set_Maximum_Iterations(const int maximum_iterations_input=0) // initially set to zero - needs to be set somewhere else!
    {maximum_iterations=maximum_iterations_input;}

    void Use_Conjugate_Gradient() // no preconditioner
    {incomplete_cholesky=false;modified_incomplete_cholesky=false;}

    void Use_Incomplete_Cholesky()
    {incomplete_cholesky=true;modified_incomplete_cholesky=false;}
    
    // use .99 for uniform grids!
    void Use_Modified_Incomplete_Cholesky(const T modified_incomplete_cholesky_coefficient_input=(T).97)
    {incomplete_cholesky=true;modified_incomplete_cholesky=true;modified_incomplete_cholesky_coefficient=modified_incomplete_cholesky_coefficient_input;} // note that both are true
    
//#####################################################################
    bool Solve(SPARSE_MATRIX_FLAT_MXN<T>& A_matrix,ARRAY<T>& x,ARRAY<T>& b,ARRAY<KRYLOV_VECTOR_BASE<T>*>& vectors,
        const T tolerance=1e-7,const bool recompute_preconditioner=true);
//#####################################################################
};
}
#endif
