//#####################################################################
// Copyright 2013.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class SCALE_MATRIX
//##################################################################### 
#ifndef __SCALE_MATRIX__
#define __SCALE_MATRIX__

#include <Core/Math_Tools/cube.h>
#include <Core/Math_Tools/pow.h>
#include <Core/Matrices/GENERAL_MATRIX_BASE.h>
#include <Core/Matrices/MATRIX.h>
#include <Core/Matrices/SYMMETRIC_MATRIX.h>
#include <Core/Matrices/ZERO_MATRIX.h>
#include <Core/Vectors/VECTOR.h>
#include <Core/Vectors/ZERO_VECTOR.h>
#include <cmath>
namespace PhysBAM{
template<class T,int d>
class SCALE_MATRIX:public GENERAL_MATRIX_BASE<T,SCALE_MATRIX<T,d>>
{
public:
    typedef T SCALAR;
    static constexpr int m=d,n=d;
    static constexpr int fixed_m=m,fixed_n=n;
    static constexpr bool is_sym=true,is_diag=true;
    static constexpr int type=mat_scale;
    T x;

    explicit SCALE_MATRIX(T z=T{}): x(z) {}

    explicit SCALE_MATRIX(INITIAL_SIZE mm,INITIAL_SIZE nn=INITIAL_SIZE(d))
        :x{}
    {}

    SCALE_MATRIX(const SCALE_MATRIX&)=default;
    SCALE_MATRIX(const IDENTITY_MATRIX<T,d>& A):x{1}{}
    SCALE_MATRIX(const ZERO_MATRIX<T,d>& A):x{}{}

    int Rows() const
    {return m;}

    int Columns() const
    {return n;}

    T operator()(const int i,const int j) const
    {assert(i==j);return x;}

    SCALE_MATRIX operator-() const
    {return SCALE_MATRIX(-x);}

    SCALE_MATRIX operator*(T a) const
    {return SCALE_MATRIX(x*a);}

    SCALE_MATRIX operator/(T a) const
    {return SCALE_MATRIX(x/a);}

    SCALE_MATRIX operator+(T a) const
    {return SCALE_MATRIX(x+a);}

    SCALE_MATRIX operator-(T a) const
    {return SCALE_MATRIX(x-a);}

    SCALE_MATRIX Transposed() const
    {return *this;}

    template<class A>
    auto Transpose_Times(const A& a) const
    {return a*x;}

    template<class A,class R=std::enable_if_t<IS_MATRIX<A>::value&&compat(n,A::fixed_n)>>
    auto Times_Transpose(const A& a) const
    {return a.Transposed()*x;}

    SCALE_MATRIX Twice_Symmetric_Part() const
    {return *this*2;}

    SCALE_MATRIX Symmetric_Part() const
    {return *this;}

    SCALE_MATRIX Diagonal_Part() const
    {return *this;}

    SCALE_MATRIX Normal_Equations_Matrix() const
    {return SCALE_MATRIX(x*x);}

    SCALE_MATRIX Outer_Product_Matrix() const
    {return SCALE_MATRIX(x*x);}

    T Trace() const
    {return x*m;}

    static SCALE_MATRIX Identity_Matrix()
    {return SCALE_MATRIX{1};}

    T Determinant() const
    {
        if constexpr(m>0) return pow<d>(x);
        else return T{};
    }

    T Max_Abs() const
    {
        if constexpr(m>0) return abs(x);
        else return T{};
    }

    template<class A>
    auto Double_Contract(const A& a)
    {return ::PhysBAM::Double_Contract(*this,a);}
};
template<int i,int j,class T,int d> auto get(const SCALE_MATRIX<T,d>& a)
{
    if constexpr(i==j) return a.x;
    return T{};
}
template<class T,int d> auto diag(const SCALE_MATRIX<T,d>& a,int i){return a.x;}
template<class T,int d> auto lower_(const SCALE_MATRIX<T,d>& a,int i,int j){return T{};}
template<class T,int d> auto upper_(const SCALE_MATRIX<T,d>& a,int i,int j){return T{};}

template<class T,int d>
inline std::ostream& operator<< (std::ostream& output_stream,const SCALE_MATRIX<T,d>& A)
{return output_stream<<(SYMMETRIC_MATRIX<T,d>()+A.x);}
}
#endif
