//#####################################################################
// Copyright 2002-2007, Ronald Fedkiw, Geoffrey Irving, Igor Neverov, Craig Schroeder, Tamar Shinar, Eftychios Sifakis, Huamin Wang, Rachel Weinstein.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class MATRIX_MXN
//#####################################################################
#ifndef __MATRIX_MXN__
#define __MATRIX_MXN__

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

template<class T>
class MATRIX_MXN:public MATRIX_BASE<T,MATRIX_MXN<T> >
{
public:
    static constexpr int fixed_m=-1,fixed_n=-1;
    typedef MATRIX_BASE<T,MATRIX_MXN<T> > BASE;
    using BASE::Times_Transpose;using BASE::Transpose_Times;
    typedef T SCALAR;using BASE::Determinant;
    int m,n; // size of the m by n matrix
    ARRAY<T> x; // one dimensional data; column major order

    MATRIX_MXN()
        :m(0),n(0)
    {}

    MATRIX_MXN(const int m_input,const int n_input);

    // Column major order
    MATRIX_MXN(const int m_input,std::initializer_list<T> list);
    MATRIX_MXN(const int m_input,const int n_input,std::initializer_list<T> list);
    explicit MATRIX_MXN(const int n_input);
    MATRIX_MXN(const INITIAL_SIZE m_input,const INITIAL_SIZE n_input);
    explicit MATRIX_MXN(const INITIAL_SIZE n_input);
    MATRIX_MXN(const MATRIX_MXN<T>& A)=default;
    template<class T2> explicit MATRIX_MXN(const MATRIX_MXN<T2>& A)
        :m(A.m),n(A.n),x(A.x,cast_data)
    {}

    template<int d>
    MATRIX_MXN(const DIAGONAL_MATRIX<T,d>& A)
        :m(d),n(d),x(m*n)
    {
        for(int i=0;i<d;i++)
            (*this)(i,i)=A(i,i);
    }
    template<int d>
    MATRIX_MXN(const SYMMETRIC_MATRIX<T,d>& A)
        :m(A.Rows()),n(A.Columns()),x(m*n)
    {
        for(int i=0;i<m;i++)
            for(int j=0;j<=i;j++)
            {
                auto a=A.Element_Lower(i,j);
                (*this)(i,j)=a;
                (*this)(j,i)=a;
            }
    }

    template<int d>
    MATRIX_MXN(const IDENTITY_MATRIX<T,d>& A)
        :m(d),n(d),x(m*n)
    {
        *this+=1;
    }

    template<int d>
    MATRIX_MXN(const SCALE_MATRIX<T,d>& A)
        :m(d),n(d),x(m*n)
    {
        *this+=A.x;
    }

    template<int mm,int nn>
        MATRIX_MXN(const ZERO_MATRIX<T,mm,nn>& A)
        :m(mm),n(nn),x(mm*nn)
    {
    }

    template<class T_MATRIX>
    explicit MATRIX_MXN(const MATRIX_BASE<T,T_MATRIX>& A)
        :m(A.Rows()),n(A.Columns()),x(m*n)
    {for(int j=0;j<n;j++) for(int i=0;i<m;i++) (*this)(i,j)=A(i,j);}

    MATRIX_MXN(MATRIX_MXN&& array)
        :m(array.m),n(array.n),x(std::move(array.x))
    {
    }

    ~MATRIX_MXN();

    int Rows() const
    {return m;}

    int Columns() const
    {return n;}

    void Resize(const int m_new,const int n_new);

    T& operator()(const int i,const int j)
    {assert((unsigned)i<(unsigned)m);assert((unsigned)j<(unsigned)n);return x(j*m+i);}

    const T& operator()(const int i,const int j) const
    {assert((unsigned)i<(unsigned)m);assert((unsigned)j<(unsigned)n);return x(j*m+i);}

    bool Valid_Index(const int i,const int j) const
    {return (unsigned)i<(unsigned)m && (unsigned)j<(unsigned)n;}

    bool operator==(const MATRIX_MXN& A) const;
    bool operator!=(const MATRIX_MXN& A) const
    {return !(*this==A);}

    MATRIX_MXN<T>& operator=(const MATRIX_MXN<T>& A);

    template<class T_MATRIX>
    MATRIX_MXN<T>& operator=(const MATRIX_BASE<T,T_MATRIX>& A)
    {if((void*)&A==(void*)this) return *this;
    x.Resize(A.Rows()*A.Columns());
    m=A.Rows();n=A.Columns();
    for(int j=0;j<n;j++) for(int i=0;i<m;i++) (*this)(i,j)=A(i,j);
    return *this;}

    template<int d> MATRIX_MXN<T>&
    operator=(const DIAGONAL_MATRIX<T,d>& A)
    {
        x.Resize(A.Rows()*A.Columns());
        m=n=d;
        x.Fill(0);
        for(int i=0;i<d;i++)
            (*this)(i,i)=A(i,i);
        return *this;
    }

    template<int d> MATRIX_MXN<T>&
    operator=(const SYMMETRIC_MATRIX<T,d>& A)
    {
        x.Resize(A.Rows()*A.Columns());
        m=n=d;
        for(int i=0;i<m;i++)
            for(int j=0;j<=i;j++)
            {
                T a=A.Element_Lower(i,j);
                (*this)(i,j)=a;
                (*this)(j,i)=a;
            }
        return *this;
    }

    MATRIX_MXN& operator=(MATRIX_MXN&& array)
    {
        m=array.m;
        n=array.n;
        x=std::move(array.x);
        return *this;
    }

    void Exchange(MATRIX_MXN& M)
    {std::swap(m,M.m);std::swap(n,M.n);x.Exchange(M.x);}

    MATRIX_MXN<T> Normal_Equations_Matrix() const;
    template<class RW> void Read(UNTYPED_ISTREAM<RW>& input);
    template<class RW> void Write(UNTYPED_OSTREAM<RW>& output) const;
//#####################################################################
};
}
#endif
