//#####################################################################
// 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.
//#####################################################################
#include <Core/Matrices/MATRIX_MXN.h>
#include <Core/Arrays/ARRAY.h>
#include <Core/Math_Tools/sqr.h>
#include <Core/Matrices/DIAGONAL_MATRIX.h>
#include <Core/Matrices/MATRIX.h>
#include <Core/Matrices/MATRIX_BASE.h>
#include <Core/Matrices/SYMMETRIC_MATRIX.h>
#include <Core/Read_Write/READ_WRITE.h>
#include <Core/Symbolics/SYMBOLIC.h>
#include <complex>
namespace PhysBAM{
//#####################################################################
// Constructor
//#####################################################################
template<class T> MATRIX_MXN<T>::
MATRIX_MXN(const int m_input,const int n_input)
    :m(m_input),n(n_input),x(m*n)
{
}
//#####################################################################
// Constructor
//#####################################################################
template<class T> MATRIX_MXN<T>::
MATRIX_MXN(const int n_input)
    :m(n_input),n(n_input),x(m*n)
{
}
//#####################################################################
// Constructor
//#####################################################################
template<class T> MATRIX_MXN<T>::
MATRIX_MXN(const INITIAL_SIZE m_input,const INITIAL_SIZE n_input)
    :m(Value(m_input)),n(Value(n_input)),x(m*n)
{
}
//#####################################################################
// Constructor
//#####################################################################
template<class T> MATRIX_MXN<T>::
MATRIX_MXN(const int m_input,std::initializer_list<T> list)
    :m(m_input),n(m_input),x(m*n)
{
    assert(list.size()==x.m);
    int k=0;
    for(const auto&t:list) x(k++)=t;
}
//#####################################################################
// Constructor
//#####################################################################
template<class T> MATRIX_MXN<T>::
MATRIX_MXN(const int m_input,const int n_input,std::initializer_list<T> list)
    :m(m_input),n(n_input),x(m*n)
{
    assert(list.size()==x.m);
    int k=0;
    for(const auto&t:list) x(k++)=t;
}
//#####################################################################
// Constructor
//#####################################################################
template<class T> MATRIX_MXN<T>::
MATRIX_MXN(const INITIAL_SIZE n_input)
    :m(Value(n_input)),n(Value(n_input)),x(m*n)
{
}
//#####################################################################
// Destructor
//#####################################################################
template<class T> MATRIX_MXN<T>::
~MATRIX_MXN()
{
}
//#####################################################################
// Function Resize
//#####################################################################
template<class T> void MATRIX_MXN<T>::
Resize(const int m_new,const int n_new)
{
    if(m_new==m && n_new==n) return;
    if(n_new==n)
    {
        x.Resize(m_new*n);
        m=m_new;
        return;
    }
    ARRAY<T> x_new(m_new*n_new);
    int m1=min(m,m_new),n1=min(n,n_new);
    for(int i=0;i<m1;i++)
        for(int j=0;j<n1;j++)
            x_new(i*n_new+j)=(*this)(i,j);
    x.Exchange(x_new);
    m=m_new;
    n=n_new;
}
//#####################################################################
// Function operator==
//#####################################################################
template<class T> bool MATRIX_MXN<T>::
operator==(const MATRIX_MXN& A) const
{
    return m==A.m && n==A.n && x==A.x;
}
//#####################################################################
// Function operator=
//#####################################################################
template<class T> MATRIX_MXN<T>& MATRIX_MXN<T>::
operator=(const MATRIX_MXN<T>& A)
{
    if(&A!=this)
    {
        x=A.x;
        m=A.m;
        n=A.n;
    }
    return *this;
}
//#####################################################################
// Function Normal_Equations_Matrix
//#####################################################################
template<class T> MATRIX_MXN<T> MATRIX_MXN<T>::
Normal_Equations_Matrix() const
{
    MATRIX_MXN<T> result(n);
    for(int j=0;j<n;j++)
        for(int i=j;i<n;i++){
            T a=0;
            for(int k=0;k<m;k++)
                a+=(*this)(k,i)*(*this)(k,j);
            result(i,j)=a;
            result(j,i)=a;
        }
    return result;
}
//#####################################################################
// Function Read
//#####################################################################
template<class T> template<class RW> void MATRIX_MXN<T>::
Read(UNTYPED_ISTREAM<RW>& input)
{
    if constexpr(is_same<RW,char>::value)
    {
        Get(input,'[');
        if(Try_Get(input,']'))
        {
            Resize(0,0);
            return;
        }
        T x;
        ARRAY<T> a;
        while(!Test(input,';') && !Test(input,']'))
        {
            Read_Binary(input,x);
            a.Append(x);
        }
        int n=a.m,m=1;
        while(Try_Get(input,';'))
        {
            m++;
            for(int i=0;i<n;i++)
            {
                Read_Binary(input,x);
                a.Append(x);
            }
        }
        Get(input,']');
        Resize(m,n);
        for(int i=0,k=0;i<m;i++)
            for(int j=0;j<n;j++)
                (*this)(i,j)=a(k++);
    }
    else Read_Binary(input,m,n,x);
}
//#####################################################################
// Function Write
//#####################################################################
template<class T> template<class RW> void MATRIX_MXN<T>::
Write(UNTYPED_OSTREAM<RW>& output) const
{
    if constexpr(is_same<RW,char>::value)
        Write_Text_Matrix(output,*this);
    else Write_Binary(output,m,n,x);
}
template class MATRIX_MXN<float>;
template class MATRIX_MXN<double>;
template class MATRIX_MXN<SYMBOLIC>;
template class MATRIX_MXN<std::complex<float>>;
template class MATRIX_MXN<std::complex<double>>;
template MATRIX_MXN<std::string >::MATRIX_MXN(int,int);
template MATRIX_MXN<std::string >::~MATRIX_MXN();
template void MATRIX_MXN<double>::Read(UNTYPED_ISTREAM<char>&);
template void MATRIX_MXN<double>::Read(UNTYPED_ISTREAM<double>&);
template void MATRIX_MXN<double>::Read(UNTYPED_ISTREAM<float>&);
template void MATRIX_MXN<double>::Write(UNTYPED_OSTREAM<char>&) const;
template void MATRIX_MXN<double>::Write(UNTYPED_OSTREAM<double>&) const;
template void MATRIX_MXN<double>::Write(UNTYPED_OSTREAM<float>&) const;
template void MATRIX_MXN<float>::Read(UNTYPED_ISTREAM<char>&);
template void MATRIX_MXN<float>::Read(UNTYPED_ISTREAM<double>&);
template void MATRIX_MXN<float>::Read(UNTYPED_ISTREAM<float>&);
template void MATRIX_MXN<float>::Write(UNTYPED_OSTREAM<char>&) const;
template void MATRIX_MXN<float>::Write(UNTYPED_OSTREAM<double>&) const;
template void MATRIX_MXN<float>::Write(UNTYPED_OSTREAM<float>&) const;
}
