//#####################################################################
// Copyright 2007-2009, Geoffrey Irving, Tamar Shinar.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Namespace MATRIX_SYSTEM
//#####################################################################
#ifndef __MATRIX_SYSTEM__
#define __MATRIX_SYSTEM__

#include <Core/Krylov_Solvers/KRYLOV_SYSTEM_BASE.h>
#include <Core/Matrices/MATRIX_FORWARD.h>
#include <Core/Utilities/TYPE_UTILITIES.h>
#include <Core/Vectors/VECTOR_FORWARD.h>
namespace PhysBAM{
//#####################################################################
// Class MATRIX_SYSTEM
//#####################################################################
template<class T_MATRIX,class T,class VECTOR_T,class T_MATRIX_PRECON=T_MATRIX>
struct MATRIX_SYSTEM:public KRYLOV_SYSTEM_BASE<T>
{
    using typename KRYLOV_SYSTEM_BASE<T>::D;
    const T_MATRIX& A;
    const T_MATRIX_PRECON* P;
    mutable VECTOR_T* temp_vector;
    ARRAY<VECTOR_T*> nullspace_vectors; // must be orthogonal and normalized
    
    MATRIX_SYSTEM(const T_MATRIX& A_input)
        :KRYLOV_SYSTEM_BASE<T>(false,true),A(A_input)
    {
    }
    virtual ~MATRIX_SYSTEM() = default;

    void Set_Preconditioner(const T_MATRIX_PRECON& preconditioner,VECTOR_T& vector)
    {
        this->use_preconditioner=true;
        P=&preconditioner;
        temp_vector=&vector;
    }
    void Multiply(const KRYLOV_VECTOR_BASE<T>& x,KRYLOV_VECTOR_BASE<T>& result,bool transpose=false) const override
    {
        const VECTOR_T& vx=dynamic_cast<const VECTOR_T&>(x);
        VECTOR_T& vresult=dynamic_cast<VECTOR_T&>(result);
        A.Times(vx.v,vresult.v);
    }
    D Inner_Product(const KRYLOV_VECTOR_BASE<T>& x,const KRYLOV_VECTOR_BASE<T>& y) const override
    {
        const VECTOR_T& vx=dynamic_cast<const VECTOR_T&>(x);
        const VECTOR_T& vy=dynamic_cast<const VECTOR_T&>(y);
        return vx.v.Dot_Double_Precision(vy.v);
    }
    void Inner_Product_Operator(KRYLOV_VECTOR_BASE<T>& x) const override {}
    T Convergence_Norm(const KRYLOV_VECTOR_BASE<T>& x) const override
    {
        const VECTOR_T& vx=dynamic_cast<const VECTOR_T&>(x);
        return vx.v.Maximum_Magnitude();
    }
    void Project(KRYLOV_VECTOR_BASE<T>& x) const override
    {
        for(auto p:nullspace_vectors)
            x.Copy(-Inner_Product(x,*p),*p,x);
    }
    void Set_Boundary_Conditions(KRYLOV_VECTOR_BASE<T>& x) const override {}
    void Project_Nullspace(KRYLOV_VECTOR_BASE<T>& x) const override {Project(x);} 
    void Apply_Preconditioner(const KRYLOV_VECTOR_BASE<T>& r,KRYLOV_VECTOR_BASE<T>& z) const override
    {
        if constexpr(is_same<SPARSE_MATRIX_FLAT_MXN<T>,T_MATRIX>::value)
        {
            const VECTOR_T& vr=dynamic_cast<const VECTOR_T&>(r);
            VECTOR_T& vz=dynamic_cast<VECTOR_T&>(z);
            P->Solve_Forward_Substitution(vr.v,temp_vector->v,true);
            P->Solve_Backward_Substitution(temp_vector->v,vz.v,false,true);
        }
        else PHYSBAM_FATAL_ERROR();
    }
};
//#####################################################################
template<class T,class TV> class KRYLOV_VECTOR_WRAPPER;
extern template struct MATRIX_SYSTEM<SPARSE_MATRIX_FLAT_MXN<float>,float,KRYLOV_VECTOR_WRAPPER<float,ARRAY<float> > >;
extern template struct MATRIX_SYSTEM<SPARSE_MATRIX_FLAT_MXN<double>,double,KRYLOV_VECTOR_WRAPPER<double,ARRAY<double> > >;
extern template struct MATRIX_SYSTEM<MATRIX_MXN<float>,float,KRYLOV_VECTOR_WRAPPER<float,ARRAY<float,int> >,MATRIX_MXN<float> >;
extern template struct MATRIX_SYSTEM<MATRIX_MXN<double>,double,KRYLOV_VECTOR_WRAPPER<double,ARRAY<double,int> >,MATRIX_MXN<double> >;
}
#endif
