//#####################################################################
// Copyright 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_BASE
//#####################################################################
#ifndef __MATRIX_BASE__
#define __MATRIX_BASE__

#include <Core/Arrays/ARRAYS_FORWARD.h>
#include <Core/Data_Structures/ELEMENT_ID.h>
#include <Core/Log/DEBUG_UTILITIES.h>
#include <Core/Math_Tools/minmax.h>
#include <Core/Math_Tools/sqr.h>
#include <Core/Matrices/GENERAL_MATRIX_BASE.h>
#include <Core/Matrices/MATRIX_FORWARD.h>
#include <Core/Matrices/MATRIX_OPERATIONS.h>
#include <Core/Read_Write/TYPED_STREAM.h>
#include <cassert>
#include <cmath>
#include <iomanip>
#include <iostream>
using std::abs;
namespace PhysBAM{

template<class T_MATRIX> struct MATRIX_INFO;
template<class T,int a,int b> struct MATRIX_INFO<MATRIX<T,a,b> >{static const int m=a,n=b;};
template<class T,int d> struct MATRIX_INFO<DIAGONAL_MATRIX<T,d> >{static const int m=d,n=d;};
template<class T,int d> struct MATRIX_INFO<SYMMETRIC_MATRIX<T,d> >{static const int m=d,n=d;};
template<class T> struct MATRIX_INFO<MATRIX_MXN<T> >{static const int m=-1,n=-1;};
template<class T> struct MATRIX_INFO<MATRIX_VIEW<T> >{static const int m=-1,n=-1;};
template<class T,class T_MATRIX> struct MATRIX_INFO<MATRIX_BASE<T,T_MATRIX> >:public MATRIX_INFO<T_MATRIX> {};

template<class T,int a,int b> struct SIZED_MATRIX_HELPER {typedef MATRIX<T,a,b> TYPE;};
template<class T,int a> struct SIZED_MATRIX_HELPER<T,a,-1> {typedef MATRIX_MXN<T> TYPE;};
template<class T,int b> struct SIZED_MATRIX_HELPER<T,-1,b> {typedef MATRIX_MXN<T> TYPE;};
template<class T> struct SIZED_MATRIX_HELPER<T,-1,-1> {typedef MATRIX_MXN<T> TYPE;};
template<class T,int a,int b> using SIZED_MATRIX = typename SIZED_MATRIX_HELPER<T,a,b>::TYPE;
template<class T_MAT> using MATRIX_COPY = SIZED_MATRIX<typename T_MAT::ELEMENT,MATRIX_INFO<T_MAT>::m,MATRIX_INFO<T_MAT>::n>;

template<class T,int a> struct SIZED_VECTOR_HELPER {typedef VECTOR<T,a> TYPE;};
template<class T> struct SIZED_VECTOR_HELPER<T,-1> {typedef ARRAY<T> TYPE;};
template<class T,int a> using SIZED_VECTOR = typename SIZED_VECTOR_HELPER<T,a>::TYPE;
template<class T_VEC> using VECTOR_COPY = SIZED_VECTOR<typename T_VEC::ELEMENT,FIXED_SIZE_VECTOR<T_VEC>::size>;

template<class T_MATRIX> struct EFFICIENT_MATRIX {static const bool value=false;};
template<class T,int m,int n> struct EFFICIENT_MATRIX<MATRIX<T,m,n> > {static const bool value=((m>=2 && m<=3 && n>=2 && n<=3) || (m==4 && n==4) || (m==0 && n==0));};
template<class T,int d> struct EFFICIENT_MATRIX<DIAGONAL_MATRIX<T,d> > {static const bool value=(d==2 || d==3);};
template<class T,int d> struct EFFICIENT_MATRIX<SYMMETRIC_MATRIX<T,d> > {static const bool value=(d==2 || d==3);};
template<class T,class T_MATRIX> struct EFFICIENT_MATRIX<MATRIX_BASE<T,T_MATRIX> >:public EFFICIENT_MATRIX<T_MATRIX>{};
template<class T,int d> struct EFFICIENT_MATRIX<VECTOR<T,d> > {static const bool value=(d<=3);};
template<class T,class T_VECTOR> struct EFFICIENT_MATRIX<ARRAY_BASE<T,T_VECTOR> >:public EFFICIENT_MATRIX<T_VECTOR>{};
template<class T> class RANDOM_NUMBERS;


namespace{
template<int line,class A,class B=void> struct ASSERT_EFFICIENT
{
    template<class T> struct EFFICIENT_OR_VOID {static const bool value=EFFICIENT_MATRIX<T>::value || is_same<T,void>::value;};
    static const bool efficient=EFFICIENT_OR_VOID<A>::value && EFFICIENT_OR_VOID<B>::value;
    struct UNUSABLE{};

    ASSERT_EFFICIENT(typename conditional<efficient,UNUSABLE,const char*>::type str)
    {}

    ASSERT_EFFICIENT(typename conditional<efficient,const char*,UNUSABLE>::type str)
    {PHYSBAM_WARNING(std::string("Base implementation used for: ")+str+".");}
};

#ifdef NDEBUG
#define WARN_IF_NOT_EFFICIENT(...)
#else
#ifdef _WIN32
#define WARN_IF_NOT_EFFICIENT(...) ASSERT_EFFICIENT<__LINE__,__VA_ARGS__>(__FUNCTION__)
#else
#define WARN_IF_NOT_EFFICIENT(...) ASSERT_EFFICIENT<__LINE__,__VA_ARGS__>(__PRETTY_FUNCTION__)
#endif
#endif
}

template<class T,class T_MATRIX>
class MATRIX_BASE:public GENERAL_MATRIX_BASE<T,T_MATRIX>
{
protected:
    MATRIX_BASE& operator=(const MATRIX_BASE&) = default;
    MATRIX_BASE& operator=(MATRIX_BASE&&) = default;
    MATRIX_BASE(const MATRIX_BASE&) = default;

public:
    static constexpr bool is_sym=false,is_diag=false;
    static constexpr int fixed_m=T_MATRIX::fixed_m;
    static constexpr int fixed_n=T_MATRIX::fixed_n;
    static constexpr int type=mat_gen;
    typedef T SCALAR;
    typedef SIZED_VECTOR<T,MATRIX_INFO<T_MATRIX>::m> LEFT_VECTOR;
    typedef SIZED_VECTOR<T,MATRIX_INFO<T_MATRIX>::n> RIGHT_VECTOR;
    typedef SIZED_VECTOR<int,MATRIX_INFO<T_MATRIX>::n> COLUMN_PERMUTATION;

    MATRIX_BASE() = default;

    ~MATRIX_BASE() = default;

    T_MATRIX& Derived()
    {return static_cast<T_MATRIX&>(*this);}

    const T_MATRIX& Derived() const
    {return static_cast<const T_MATRIX&>(*this);}

    T operator()(int i,int j) const
    {return Derived()(i,j);}

    T& operator()(int i,int j)
    {return Derived()(i,j);}

    int Rows() const
    {return Derived().Rows();}

    int Columns() const
    {return Derived().Columns();}

protected:
    // Logic only; no type or bounds checking; requires all arguments to be distinct
    template<class T_MATRIX0,class T_MATRIX1,class T_MATRIX2>
    static void Add_Times_Matrix_Helper(const T_MATRIX0& A,const T_MATRIX1& B,T_MATRIX2& C)
    {for(int j=0;j<B.Columns();j++) for(int k=0;k<A.Columns();k++) for(int i=0;i<A.Rows();i++) C(i,j)+=A(i,k)*B(k,j);}

    template<class T_MATRIX0,class T_MATRIX1,class T_MATRIX2>
    static void Add_Transpose_Times_Matrix_Helper(const T_MATRIX0& A,const T_MATRIX1& B,T_MATRIX2& C)
    {for(int j=0;j<B.Columns();j++) for(int i=0;i<A.Columns();i++) for(int k=0;k<A.Rows();k++) C(i,j)+=A(k,i)*B(k,j);}

    template<class T_MATRIX0,class T_MATRIX1,class T_MATRIX2>
    static void Add_Times_Transpose_Matrix_Helper(const T_MATRIX0& A,const T_MATRIX1& B,T_MATRIX2& C)
    {for(int k=0;k<A.Columns();k++) for(int j=0;j<B.Rows();j++) for(int i=0;i<A.Rows();i++) C(i,j)+=A(i,k)*B(j,k);}

    template<class T_MATRIX0,class T_VECTOR,class T_VECTOR1>
    static void Add_Transpose_Times_Vector_Helper(const T_MATRIX0& A,const ARRAY_BASE<T,T_VECTOR>& v,ARRAY_BASE<T,T_VECTOR1>& y)
    {for(int i=0;i<y.Size();i++) for(int j=0;j<v.Size();j++) y(i)+=A(j,i)*v(j);}

    template<class T_MATRIX0,class T_VECTOR,class T_VECTOR1>
    static void Add_Times_Vector_Helper(const T_MATRIX0& A,const ARRAY_BASE<T,T_VECTOR>& v,ARRAY_BASE<T,T_VECTOR1>& y)
    {for(int j=0;j<v.Size();j++) for(int i=0;i<y.Size();i++) y(i)+=A(i,j)*v(j);}

    template<class T_MATRIX0,class T_VECTOR,class T_VECTOR1>
    static void Subtract_Times_Vector_Helper(const T_MATRIX0& A,const ARRAY_BASE<T,T_VECTOR>& v,ARRAY_BASE<T,T_VECTOR1>& y)
    {for(int j=0;j<v.Size();j++) for(int i=0;i<y.Size();i++) y(i)-=A(i,j)*v(j);}

    template<class T_MATRIX0,class T_MATRIX1>
    static bool Test_Aliased(const T_MATRIX0& A,const T_MATRIX1& B)
    {assert((void*)&A!=(void*)&B);return false;}

    template<class T_MATRIX0>
    static bool Test_Aliased(const T_MATRIX0& A,const T_MATRIX0& B)
    {return &A==&B;}
public:

    // With bounds checking; arguments may be the same.
    template<class T_MATRIX0,class T_MATRIX1,class T_MATRIX2>
    static void Add_Times(const T_MATRIX0& A,const T_MATRIX1& B,MATRIX_BASE<T,T_MATRIX2>& C)
    {WARN_IF_NOT_EFFICIENT(T_MATRIX0,T_MATRIX1);assert(A.Columns()==B.Rows() && C.Rows()==A.Rows() && C.Columns()==B.Columns());
    if(Test_Aliased(A,C) || Test_Aliased(B,C)){T_MATRIX2 M(INITIAL_SIZE(C.Rows()),INITIAL_SIZE(C.Columns()));
    Add_Times_Matrix_Helper(A,B,M);C.Derived()+=M;}else Add_Times_Matrix_Helper(A,B,C);}

    template<class T_MATRIX0,class T_MATRIX1,class T_MATRIX2>
    static void Add_Transpose_Times(const T_MATRIX0& A,const T_MATRIX1& B,MATRIX_BASE<T,T_MATRIX2>& C)
    {WARN_IF_NOT_EFFICIENT(T_MATRIX0,T_MATRIX1);assert(A.Rows()==B.Rows() && C.Rows()==A.Columns() && C.Columns()==B.Columns());
    if(Test_Aliased(A,C) || Test_Aliased(B,C)){T_MATRIX2 M(INITIAL_SIZE(C.Rows()),INITIAL_SIZE(C.Columns()));
    Add_Transpose_Times_Matrix_Helper(A,B,M);C.Derived()+=M;}else Add_Transpose_Times_Matrix_Helper(A,B,C);}

    template<class T_MATRIX0,class T_MATRIX1,class T_MATRIX2>
    static void Add_Times_Transpose(const T_MATRIX0& A,const T_MATRIX1& B,MATRIX_BASE<T,T_MATRIX2>& C)
    {WARN_IF_NOT_EFFICIENT(T_MATRIX0,T_MATRIX1);assert(A.Columns()==B.Columns() && C.Rows()==A.Rows() && C.Columns()==B.Rows());
    if(Test_Aliased(A,C) || Test_Aliased(B,C)){T_MATRIX2 M(INITIAL_SIZE(C.Rows()),INITIAL_SIZE(C.Columns()));
    Add_Times_Transpose_Matrix_Helper(A,B,M);C.Derived()+=M;}else Add_Times_Transpose_Matrix_Helper(A,B,C);}

    template<class T_MATRIX0,class T_VECTOR,class T_VECTOR1>
    static void Add_Transpose_Times(const MATRIX_BASE<T,T_MATRIX0>& A,const ARRAY_BASE<T,T_VECTOR>& v,ARRAY_BASE<T,T_VECTOR1>& y)
    {WARN_IF_NOT_EFFICIENT(T_MATRIX0,T_VECTOR1);assert(A.Rows()==v.Size() && A.Columns()==y.Size());
    if(Test_Aliased(v,y)){VECTOR_COPY<T_VECTOR> u(v);Add_Transpose_Times_Vector_Helper(A,u,y);}else Add_Transpose_Times_Vector_Helper(A,v,y);}

    template<class T_MATRIX0,class T_VECTOR,class T_VECTOR1>
    static void Add_Times(const MATRIX_BASE<T,T_MATRIX0>& A,const ARRAY_BASE<T,T_VECTOR>& v,ARRAY_BASE<T,T_VECTOR1>& y)
    {WARN_IF_NOT_EFFICIENT(T_MATRIX0,T_VECTOR);assert(A.Columns()==v.Size() && A.Rows()==y.Size());
    if(Test_Aliased(v,y)){VECTOR_COPY<T_VECTOR> u(v);Add_Times_Vector_Helper(A,u,y);}else Add_Times_Vector_Helper(A,v,y);}

    template<class T_MATRIX0,class T_VECTOR,class T_VECTOR1>
    static void Subtract_Times(const MATRIX_BASE<T,T_MATRIX0>& A,const ARRAY_BASE<T,T_VECTOR>& v,ARRAY_BASE<T,T_VECTOR1>& y)
    {WARN_IF_NOT_EFFICIENT(T_MATRIX0,T_VECTOR);assert(A.Columns()==v.Size() && A.Rows()==y.Size());
    if(Test_Aliased(v,y)){VECTOR_COPY<T_VECTOR> u(v);Subtract_Times_Vector_Helper(A,u,y);}else Subtract_Times_Vector_Helper(A,v,y);}

    template<class T_VECTOR,class T_VECTOR1>
    void Times(const ARRAY_BASE<T,T_VECTOR>& y,ARRAY_BASE<T,T_VECTOR1>& z) const
    {z.Fill(T());Add_Times(Derived(),y.Derived(),z.Derived());}

    T Infinity_Norm() const
    {
        T max_sum=0;
        for(int i=0;i<Rows();i++)
        {
            T sum=0;
            for(int j=0;j<Columns();j++)
                sum+=abs((*this)(i,j));
            max_sum=max(sum,max_sum);
        }
        return max_sum;
    }

    template<class T_MATRIX1>
    void Set_Submatrix(int istart,int jstart,const MATRIX_BASE<T,T_MATRIX1>& a)
    {for(int i=0;i<a.Rows();i++) for(int j=0;j<a.Columns();j++) (*this)(istart+i,jstart+j)=a(i,j);}

    template<class MAT> std::enable_if_t<(IS_MATRIX<MAT>::value&&MAT::fixed_m>=0)>
    Set_Submatrix(int istart,int jstart,const MAT& a)
    {
        []<int... i>(auto& t,const auto& a,int j,int k,std::integer_sequence<int,i...>)
        {
            ((t(j+i%MAT::m,k+i/MAT::m)=get<i%MAT::m,i/MAT::m>(a)),...);
        }(*this,a,istart,jstart,std::make_integer_sequence<int,MAT::m*MAT::n>{});
    }

    template<class T_MATRIX1>
    void Add_To_Submatrix(int istart,int jstart,const MATRIX_BASE<T,T_MATRIX1>& a)
    {for(int i=0;i<a.Rows();i++) for(int j=0;j<a.Columns();j++) (*this)(istart+i,jstart+j)+=a(i,j);}

    template<class MAT> std::enable_if_t<(IS_MATRIX<MAT>::value&&MAT::fixed_m>=0)>
    Add_To_Submatrix(int istart,int jstart,const MAT& a)
    {
        []<int... i>(auto& t,const auto& a,int j,int k,std::integer_sequence<int,i...>)
        {
            ((t(j+i%MAT::m,k+i/MAT::m)+=get<i%MAT::m,i/MAT::m>(a)),...);
        }(*this,a,istart,jstart,std::make_integer_sequence<int,MAT::m*MAT::n>{});
    }

    template<int d>
    void Add_To_Submatrix(int istart,int jstart,const DIAGONAL_MATRIX<T,d>& a)
    {
        []<int... i>(auto& t,const auto& a,int j,int k,std::integer_sequence<int,i...>)
        {
            ((t(j+i,k+i)+=a(i,i)),...);
        }(*this,a,istart,jstart,std::make_integer_sequence<int,d>{});
    }

    template<int mm,int nn>
    void Get_Submatrix(int istart,int jstart,MATRIX<T,mm,nn>& a) const
    {
        []<int... i>(const auto& t,auto& a,int j,int k,std::integer_sequence<int,i...>)
        {
            ((a(i%mm,i/mm)=t(j+i%mm,k+i/mm)),...);
        }(*this,a,istart,jstart,std::make_integer_sequence<int,mm*nn>{});
    }

    template<class T_MATRIX1>
    void Get_Submatrix(int istart,int jstart,MATRIX_BASE<T,T_MATRIX1>& a) const
    {for(int i=0;i<a.Rows();i++) for(int j=0;j<a.Columns();j++) a(i,j)=(*this)(istart+i,jstart+j);}

    auto Get_Column(int j) const
    {return ::PhysBAM::Get_Column(Derived(),j);}

    auto Get_Row(int i) const
    {return ::PhysBAM::Get_Row(Derived(),i);}

    auto Column(int j) const
    {return ::PhysBAM::Get_Column(Derived(),j);}

    auto Row(int i) const
    {return ::PhysBAM::Get_Row(Derived(),i);}

    template<class T_VECTOR>
    void Set_Column(int j,const ARRAY_BASE<T,T_VECTOR>& a)
    {::PhysBAM::Set_Column(Derived(),j,a);}

    template<class T_VECTOR>
    void Set_Row(int i,const ARRAY_BASE<T,T_VECTOR>& a)
    {::PhysBAM::Set_Row(Derived(),i,a);}

    template<class T_VECTOR>
    void Add_Column(int j,const ARRAY_BASE<T,T_VECTOR>& a)
    {::PhysBAM::Add_Column(Derived(),j,a);}

    template<class T_VECTOR>
    void Add_Row(int i,const ARRAY_BASE<T,T_VECTOR>& a)
    {::PhysBAM::Add_Row(Derived(),i,a);}

    void Set_Zero_Matrix()
    {for(int i=0;i<Rows();i++) for(int j=0;j<Columns();j++) (*this)(i,j)=T();}

    void Set_Identity_Matrix()
    {Set_Zero_Matrix();Derived()+=1;}

    static T_MATRIX Identity_Matrix(int n)
    {T_MATRIX A(n);A+=1;return A;}

    template<class T_VECTOR>
    RIGHT_VECTOR Lower_Triangular_Solve(const ARRAY_BASE<T,T_VECTOR>& b) const
    {assert(Rows()==Columns() && Columns()==b.Size());RIGHT_VECTOR x(INITIAL_SIZE(b.Size()));
    for(int i=0;i<Columns();i++){
        x(i)=b(i);
        for(int j=0;j<i;j++) x(i)-=(*this)(i,j)*x(j);
        x(i)/=(*this)(i,i);}
    return x;}

    template<class T_VECTOR>
    LEFT_VECTOR Transpose_Lower_Triangular_Solve(const ARRAY_BASE<T,T_VECTOR>& b) const
    {assert(Rows()==Columns() && Columns()==b.Size());LEFT_VECTOR x(b);
    for(int i=0;i<Columns();i++){
        x(i)/=(*this)(i,i);
        for(int j=i+1;j<Columns();j++) x(j)-=(*this)(i,j)*x(i);}
    return x;}

    template<class T_VECTOR>
    RIGHT_VECTOR Upper_Triangular_Solve(const ARRAY_BASE<T,T_VECTOR>& b) const
    {assert(Rows()==Columns() && Columns()==b.Size());RIGHT_VECTOR x(INITIAL_SIZE(b.Size()));
    for(int i=Columns()-1;i>=0;i--){
        x(i)=b(i);
        for(int j=Columns()-1;j>i;j--) x(i)-=(*this)(i,j)*x(j);
        x(i)/=(*this)(i,i);}
    return x;}

    template<class T_VECTOR>
    LEFT_VECTOR Transpose_Upper_Triangular_Solve(const ARRAY_BASE<T,T_VECTOR>& b) const
    {assert(Rows()==Columns() && Columns()==b.Size());LEFT_VECTOR x(b);
    for(int i=Columns()-1;i>=0;i--){
        x(i)/=(*this)(i,i);
        for(int j=i-1;j>=0;j--) x(j)-=(*this)(i,j)*x(i);}
    return x;}

    template<class T_MATRIX1>
    T_MATRIX1 Upper_Triangular_Solve(const MATRIX_BASE<T,T_MATRIX1>& b) const
    {assert(Rows()==Columns() && Columns()==b.Rows());T_MATRIX1 x(INITIAL_SIZE(b.Rows()),INITIAL_SIZE(b.Columns()));
    for(int bcol=0;bcol<b.Columns();bcol++)
        for(int i=Columns()-1;i>=0;i--){
            x(i,bcol)=b(i,bcol);
            for(int j=Columns()-1;j>=i+1;j--) x(i,bcol)-=(*this)(i,j)*x(j,bcol);
            x(i,bcol)/=(*this)(i,i);}
    return x;}

    template<class T_VECTOR>
    RIGHT_VECTOR In_Place_Cholesky_Solve(const ARRAY_BASE<T,T_VECTOR>& b)
    {assert(Rows()==Columns());In_Place_Cholesky_Factorization();return Transpose_Upper_Triangular_Solve(Lower_Triangular_Solve(b));}

    template<class T_VECTOR>
    RIGHT_VECTOR Cholesky_Solve(const ARRAY_BASE<T,T_VECTOR>& b) const
    {return T_MATRIX(Derived()).In_Place_Cholesky_Solve(b);}

    void In_Place_Cholesky_Inverse()
    {return T_MATRIX(Derived()).In_Place_Cholesky_Inverse(*this);}

    template<class T_MATRIX1>
    void Cholesky_Inverse(MATRIX_BASE<T,T_MATRIX1>& inverse) const
    {return T_MATRIX(Derived()).In_Place_Cholesky_Inverse(inverse);}

    template<class T_VECTOR,class T_MATRIX1>
    static typename T_MATRIX1::LEFT_VECTOR Householder_Transform(const ARRAY_BASE<T,T_VECTOR>& b,const MATRIX_BASE<T,T_MATRIX1>& V)
    {assert(V.Rows()==b.Size());typename T_MATRIX1::LEFT_VECTOR result(b),v((INITIAL_SIZE)V.Rows());
    for(int j=0;j<V.Columns();j++){result=result.Householder_Transform(V.Get_Column(j));}
    return result;}

    template<class T_VECTOR>
    LEFT_VECTOR Householder_QR_Solve(const ARRAY_BASE<T,T_VECTOR>& b) const
    {T_MATRIX V,R;Householder_QR_Factorization(V,R);LEFT_VECTOR c=Householder_Transform(b,V),c_short((INITIAL_SIZE)Columns());
    for(int i=0;i<Columns();i++) c_short(i)=c(i);
    return R.Upper_Triangular_Solve(c_short);}

    template<class T_VECTOR>
    RIGHT_VECTOR In_Place_PLU_Solve(const ARRAY_BASE<T,T_VECTOR>& b)
    {assert(Rows()==Columns());COLUMN_PERMUTATION p;T_MATRIX L;In_Place_PLU_Factorization(L,p);
    return Upper_Triangular_Solve(L.Lower_Triangular_Solve(b.Subset(p)));}

    template<class T_VECTOR>
    RIGHT_VECTOR PLU_Solve(const ARRAY_BASE<T,T_VECTOR>& b) const
    {return T_MATRIX(Derived()).In_Place_PLU_Solve(b);}

    template<class T_MATRIX1>
    void PLU_Inverse(MATRIX_BASE<T,T_MATRIX1>& inverse) const
    {(T_MATRIX1(Derived())).In_Place_PLU_Inverse(inverse);}

    template<class T_VECTOR>
    RIGHT_VECTOR In_Place_LU_Solve(const ARRAY_BASE<T,T_VECTOR>& b)
    {assert(Rows()==Columns());T_MATRIX L;In_Place_LU_Factorization(L);return Upper_Triangular_Solve(L.Lower_Triangular_Solve(b));}

    template<class T_VECTOR>
    RIGHT_VECTOR LU_Solve(const ARRAY_BASE<T,T_VECTOR>& b) const
    {return T_MATRIX(Derived()).In_Place_LU_Solve(b);}

    template<class T_MATRIX1>
    void LU_Inverse(MATRIX_BASE<T,T_MATRIX1>& inverse) const
    {(T_MATRIX1(Derived())).In_Place_LU_Inverse(inverse);}

    template<class T_VECTOR>
    RIGHT_VECTOR Inverse_Times(const ARRAY_BASE<T,T_VECTOR>& b) const
    {return PLU_Solve(b);}

//#####################################################################
    template<class T_MATRIX1> void In_Place_Cholesky_Inverse(MATRIX_BASE<T,T_MATRIX1>& inverse);
    template<class T_MATRIX1> void In_Place_PLU_Inverse(MATRIX_BASE<T,T_MATRIX1>& inverse);
    template<class T_MATRIX1> void In_Place_LU_Inverse(MATRIX_BASE<T,T_MATRIX1>& inverse);
    template<class T_MATRIX1,class T_MATRIX2> void Householder_QR_Factorization(MATRIX_BASE<T,T_MATRIX1>& V,MATRIX_BASE<T,T_MATRIX2>& R) const;
    template<class T_VECTOR0,class T_VECTOR1> void In_Place_Robust_Householder_QR_Solve(ARRAY_BASE<T,T_VECTOR0>& b,ARRAY_BASE<int,T_VECTOR1>& p); // this=Q
    template<class T_MATRIX1> void In_Place_PLU_Factorization(MATRIX_BASE<T,T_MATRIX1>& L,COLUMN_PERMUTATION& p); // this=U
    void In_Place_Cholesky_Factorization(); // this=L
    template<class T_MATRIX1> void In_Place_LU_Factorization(MATRIX_BASE<T,T_MATRIX1>& L); // this=U
    void Fill_Uniform(RANDOM_NUMBERS<T>& rand,T a,T b);
protected:
    T Determinant() const;
//#####################################################################
};

template<int i,int j,class T,class T_MATRIX> auto get(const MATRIX_BASE<T,T_MATRIX>& a){return a(i,j);}
template<class T,class T_MATRIX> auto get(const MATRIX_BASE<T,T_MATRIX>& a,int i,int j){return a(i,j);}
template<class T,class T_MATRIX> auto diag(const MATRIX_BASE<T,T_MATRIX>& a,int i){return a(i,i);}
template<class T,class T_MATRIX> auto lower(const MATRIX_BASE<T,T_MATRIX>& a,int i,int j){return a(i,j);}
template<class T,class T_MATRIX> auto lower_(const MATRIX_BASE<T,T_MATRIX>& a,int i,int j){return a(i,j);}
template<class T,class T_MATRIX> auto upper(const MATRIX_BASE<T,T_MATRIX>& a,int i,int j){return a(i,j);}
template<class T,class T_MATRIX> auto upper_(const MATRIX_BASE<T,T_MATRIX>& a,int i,int j){return a(i,j);}

template<class T,class T_MATRIX> inline std::ostream& operator<<(std::ostream& out,const MATRIX_BASE<T,T_MATRIX>& A)
{
    RAW_OSTREAM r(out);
    UNTYPED_OSTREAM<char> s(r);
    A.Derived().Write(s);
    return out;
}
template<class T,class T_MATRIX> inline std::istream& operator>>(std::istream& in,MATRIX_BASE<T,T_MATRIX>& A)
{
    RAW_ISTREAM r(in);
    UNTYPED_ISTREAM<char> s(r);
    A.Derived().Read(s);
    return in;
}

template<class T,class T_MATRIX>
T_MATRIX exp(const NO_MATRIX& M);
template<class T_MATRIX> std::enable_if_t<IS_MATRIX<T_MATRIX>::value>
Write_Text_Matrix(UNTYPED_OSTREAM<char>& output,const T_MATRIX& m);
template<class T_MATRIX> std::enable_if_t<IS_MATRIX<T_MATRIX>::value>
Write_Text_Matrix_Upper(UNTYPED_OSTREAM<char>& output,const T_MATRIX& m);
template<class T_MATRIX> std::enable_if_t<IS_MATRIX<T_MATRIX>::value>
Read_Text_Matrix(UNTYPED_ISTREAM<char>& input,T_MATRIX& m);
template<class T_MATRIX> std::enable_if_t<IS_MATRIX<T_MATRIX>::value>
Read_Text_Matrix_Upper(UNTYPED_ISTREAM<char>& input,T_MATRIX& m);
template<class RW,class T,class T_MATRIX,int s> std::enable_if_t<IS_MATRIX<T_MATRIX>::value>
Write_Fixed_Matrix(UNTYPED_OSTREAM<RW>& output,const T_MATRIX& m,T(& x)[s]);
template<class RW,class T,class T_MATRIX,int s> std::enable_if_t<IS_MATRIX<T_MATRIX>::value>
Write_Fixed_Matrix_Upper(UNTYPED_OSTREAM<RW>& output,const T_MATRIX& m,T(& x)[s]);
template<class RW,class T,class T_MATRIX,int s> std::enable_if_t<IS_MATRIX<T_MATRIX>::value>
Read_Fixed_Matrix(UNTYPED_ISTREAM<RW>& input,T_MATRIX& m,T(& x)[s]);
template<class RW,class T,class T_MATRIX,int s> std::enable_if_t<IS_MATRIX<T_MATRIX>::value>
Read_Fixed_Matrix_Upper(UNTYPED_ISTREAM<RW>& input,T_MATRIX& m,T(& x)[s]);
}
#endif
