//#####################################################################
// Copyright 2007, Jon Gretarsson, Avi Robinson-Mosher, Tamar Shinar, Jonathan Su.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class SPARSE_MATRIX_FLAT_MXN
//#####################################################################
#ifndef __SPARSE_MATRIX_FLAT_MXN__
#define __SPARSE_MATRIX_FLAT_MXN__

#include <Core/Arrays/ARRAY.h>
namespace PhysBAM{

template<class T> class INTERVAL;

template<class T>
struct SPARSE_MATRIX_ENTRY
{
    int j;T a;
    SPARSE_MATRIX_ENTRY():j(-1),a(0){}
    SPARSE_MATRIX_ENTRY(int index,T value):j(index),a(value){}
    bool operator<(const SPARSE_MATRIX_ENTRY& s) const {return j<s.j;}

    template<class RW> void Read(UNTYPED_ISTREAM<RW>& input)
    {Read_Binary_Group(input,j,a);}

    template<class RW> void Write(UNTYPED_OSTREAM<RW>& output) const
    {Write_Binary_Group(output,j,a);}
};

// Assumes sorted rows:
// operator(), Find_Index

template<class T>
class SPARSE_MATRIX_FLAT_MXN
{
public:
    typedef T SCALAR;
    int m,n; // size of the m by n matrix
    ARRAY<int> offsets;
    ARRAY<SPARSE_MATRIX_ENTRY<T> > A;
    SPARSE_MATRIX_FLAT_MXN<T>* Q;
    SPARSE_MATRIX_FLAT_MXN<T>* L;
    SPARSE_MATRIX_FLAT_MXN<T>* C;
    ARRAY<int> diagonal_index;

    SPARSE_MATRIX_FLAT_MXN();

    SPARSE_MATRIX_FLAT_MXN(const SPARSE_MATRIX_FLAT_MXN<T>& matrix);

    ~SPARSE_MATRIX_FLAT_MXN();

    SPARSE_MATRIX_FLAT_MXN& operator=(const SPARSE_MATRIX_FLAT_MXN& matrix)
    {m=matrix.m;n=matrix.n;offsets=matrix.offsets;A=matrix.A;diagonal_index=matrix.diagonal_index;delete C;C=0;
    return *this;}

    T operator()(int i,int j) const
    {int index=Find_Index(i,j);assert(A(index).j==j);return A(index).a;}

    T& operator()(int i,int j);

    // void f(int r,int c,T a);
    template<class F>
    void For_Each(F f) const
    {
        for(int i=0,b=offsets(0);i<m;i++)
        {
            int e=offsets(i+1);
            for(int j=b;j<e;j++)
                f(i,A(j).j,A(j).a);
            b=e;
        }
    }

    // void f(int r,int c,T& a);
    template<class F>
    void For_Each(F f)
    {
        for(int i=0,b=offsets(0);i<m;i++)
        {
            int e=offsets(i+1);
            for(int j=b;j<e;j++)
                f(i,A(j).j,A(j).a);
            b=e;
        }
    }

    template<class B,class F,class E>
    void For_Each(B begin_row,F f,E end_row) const
    {
        for(int i=0,b=offsets(0);i<m;i++)
        {
            int e=offsets(i+1);
            begin_row(i);
            for(int j=b;j<e;j++)
                f(i,A(j).j,A(j).a);
            end_row(i);
            b=e;
        }
    }
    
//#####################################################################
    SPARSE_MATRIX_FLAT_MXN<T>* Create_Submatrix(const INTERVAL<int>& rows);
    void Set_Row_Lengths(ARRAY_VIEW<int> lengths);
    int Find_Index(int i,int j) const;
    bool Element_Present(int i,int j) const;
    void Times(ARRAY_VIEW<const T> x,ARRAY_VIEW<T> result) const;
    void Transpose_Times(ARRAY_VIEW<const T> x,ARRAY_VIEW<T> result) const;
    void Times_Add(ARRAY_VIEW<const T> x,ARRAY_VIEW<T> result) const;
    void Times_Add_Row(ARRAY_VIEW<const T> x,ARRAY_VIEW<T> result,int row) const;
    void Times_Subtract(ARRAY_VIEW<const T> x,ARRAY_VIEW<T> result) const;
    void Times_Threaded(ARRAY_VIEW<const T> x,ARRAY_VIEW<T> result) const;
    void Times_Add_Threaded(ARRAY_VIEW<const T> x,ARRAY_VIEW<T> result) const;
    void Times_Subtract_Threaded(ARRAY_VIEW<const T> x,ARRAY_VIEW<T> result) const;
    void Transpose_Times_Add(ARRAY_VIEW<const T> x,ARRAY_VIEW<T> result) const;
    void Transpose_Times_Subtract(ARRAY_VIEW<const T> x,ARRAY_VIEW<T> result) const;
    void Negate();
    SPARSE_MATRIX_FLAT_MXN<T>& operator*=(T a);
    SPARSE_MATRIX_FLAT_MXN<T>& operator/=(T a);
    SPARSE_MATRIX_FLAT_MXN<T>& operator+=(T a);
    SPARSE_MATRIX_FLAT_MXN<T>& operator-=(T a);
    void Compress(SPARSE_MATRIX_FLAT_MXN<T>& compressed) const;
    void Transpose(SPARSE_MATRIX_FLAT_MXN<T>& A_transpose) const;
    SPARSE_MATRIX_FLAT_MXN<T> Times_Transpose(const SPARSE_MATRIX_FLAT_MXN<T>& rhs);
    SPARSE_MATRIX_FLAT_MXN<T> Times_Diagonal_Times(ARRAY_VIEW<const T> diagonal,const SPARSE_MATRIX_FLAT_MXN<T>& rhs); // (*this) * diagonal * (rhs)
    SPARSE_MATRIX_FLAT_MXN<T> Scale_Rows(ARRAY_VIEW<const T> d) const;
    SPARSE_MATRIX_FLAT_MXN<T> operator+(const SPARSE_MATRIX_FLAT_MXN<T>& A_rhs) const;
    SPARSE_MATRIX_FLAT_MXN<T> operator-(const SPARSE_MATRIX_FLAT_MXN<T>& A_rhs) const;
    SPARSE_MATRIX_FLAT_MXN<T> operator*(const SPARSE_MATRIX_FLAT_MXN<T>& rhs) const;
    void Set_Times_Diagonal(ARRAY_VIEW<const T> D);
    void Set_Diagonal_Times(ARRAY_VIEW<const T> D);
    void Reset(int c);
    void Append_Entry_To_Current_Row(int c,T a);
    void Finish_Row();
    void Sort_Entries();
    void Get_Row(ARRAY<SPARSE_MATRIX_ENTRY<T> >& q_i,int row);
    void Fast_Sparse_Multiply(ARRAY<SPARSE_MATRIX_ENTRY<T> >& q,ARRAY<SPARSE_MATRIX_ENTRY<T> >& l);
    void Initialize_Diagonal_Index();
    void Solve_Forward_Substitution(ARRAY_VIEW<const T> b,ARRAY_VIEW<T> x,const bool diagonal_is_identity=false,const bool diagonal_is_inverted=false) const;
    void Solve_Backward_Substitution(ARRAY_VIEW<const T> b,ARRAY_VIEW<T> x,const bool diagonal_is_identity=false,const bool diagonal_is_inverted=false) const;
    // actually an LU saving square roots, with an inverted diagonal saving divides
    void Construct_Incomplete_Cholesky_Factorization(const bool modified_version=true,T modified_coefficient=.97,T zero_tolerance=1e-8,T zero_replacement=1e-8);
    // actually an LU saving square roots, with an inverted diagonal saving divides
    void In_Place_Incomplete_Cholesky_Factorization(const bool modified_version=true,T modified_coefficient=.97,T zero_tolerance=1e-8,T zero_replacement=1e-8);
    void Gauss_Seidel_Single_Iteration(ARRAY_VIEW<T> x,ARRAY_VIEW<const T> b);
    void Gauss_Seidel_Solve(ARRAY_VIEW<T> x,ARRAY_VIEW<const T> b,T tolerance=1e-12,int max_iterations=1000000);
    bool Positive_Diagonal_And_Nonnegative_Row_Sum(T tolerance=1e-7) const;
    void Conjugate_With_Diagonal_Matrix(ARRAY_VIEW<T> x);

    // Does not remove rows from the matrix; just removes zero entries from A.
    void Prune_Zero_Entries();
    template<class RW> void Read(UNTYPED_ISTREAM<RW>& input);
    template<class RW> void Write(UNTYPED_OSTREAM<RW>& output);
//#####################################################################
};
template<class T> std::ostream& operator<<(std::ostream& output_stream,const SPARSE_MATRIX_FLAT_MXN<T>& A);
}
#endif
