//#####################################################################
// Copyright 2019, Craig Schroeder.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.text.
//#####################################################################
#include <Core/Matrices/LAPACK.h>
namespace PhysBAM{

// Destroys M
template<class T>
void Linear_Solve_Raw_(MATRIX_VIEW<T> M,
    MATRIX_VIEW<const T> rhs,MATRIX_VIEW<T> sol)
{
#ifdef USE_LAPACK
    assert(M.Rows()==rhs.Rows());
    assert(M.Columns()==sol.Rows());
    assert(M.Rows()==M.Columns());
    assert(rhs.Columns()==sol.Columns());
    sol=rhs;

    ARRAY<int> jpvt(M.n);
    int rnk = 0;

    int ret = 0;
    T rcond=std::numeric_limits<T>::epsilon();
    if constexpr(is_same_v<T,float>)
        ret=LAPACKE_sgesv(LAPACK_COL_MAJOR,M.m,rhs.Columns(),M.x,M.s,
            jpvt.Get_Array_Pointer(),sol.x,sol.s,rcond,&rnk);
    else if constexpr(is_same_v<T,double>)
        ret=LAPACKE_dgesv(LAPACK_COL_MAJOR,M.m,rhs.Columns(),M.x,M.s,
            jpvt.Get_Array_Pointer(),sol.x,sol.s,rcond,&rnk);
    else PHYSBAM_FATAL_ERROR("LAPACK only defined for float/double types");

    PHYSBAM_ASSERT(!ret);
#else
    PHYSBAM_FATAL_ERROR("This requies USE_LAPACK");
#endif
}


// Destroys M
template<class T>
void Least_Squares_Solve_Raw_(MATRIX_VIEW<T> M,
    MATRIX_VIEW<const T> rhs,MATRIX_VIEW<T> sol)
{
#ifdef USE_LAPACK
    assert(M.Rows()==rhs.Rows());
    assert(M.Columns()==sol.Rows());
    assert(rhs.Columns()==sol.Columns());

    MATRIX_MXN<T> TS;
    MATRIX_VIEW<T> VS(sol.Derived());

    bool b=false;
    if(sol.Rows()<rhs.Rows())
    {
        TS.Resize(rhs.Rows(),rhs.Columns());
        VS.Set(TS);
        b=true;
    }
    VS.m=rhs.Rows();
    VS=rhs;
    VS.m=sol.Rows();

    ARRAY<int> jpvt(M.n);
    int rnk = 0;
    int ret = 0;
    T rcond=std::numeric_limits<T>::epsilon();

    if constexpr(is_same_v<T,float>)
        ret=LAPACKE_sgelsy(LAPACK_COL_MAJOR,M.m,M.n,rhs.Columns(),M.x,M.s,
            VS.x,VS.s,jpvt.Get_Array_Pointer(),rcond,&rnk);
    else if constexpr(is_same_v<T,double>)
        ret=LAPACKE_dgelsy(LAPACK_COL_MAJOR,M.m,M.n,rhs.Columns(),M.x,M.s,
            VS.x,VS.s,jpvt.Get_Array_Pointer(),rcond,&rnk);
    else PHYSBAM_FATAL_ERROR("LAPACK only defined for float/double types");

    PHYSBAM_ASSERT(!ret);
    if(b) sol.Derived()=VS;
#else
    PHYSBAM_FATAL_ERROR("This requies USE_LAPACK");
#endif
}

// Compute (complex) eigenvalues for non-symmetric real-valued matrix.
// Destroys M
template<class T>
void Compute_Eigenvalues_Raw_(MATRIX_VIEW<T> M,ARRAY_VIEW<std::complex<T>> eigenvalues)
{
#ifdef USE_LAPACK
    PHYSBAM_ASSERT(M.m==M.n);
    ARRAY<T> r(eigenvalues.m),i(eigenvalues.m);

    int ret=0;
    if constexpr(is_same_v<T,float>)
        ret=LAPACKE_sgeev(LAPACK_COL_MAJOR,'N','N',M.m,M.x,M.s,
            r.Get_Array_Pointer(),i.Get_Array_Pointer(),0,1,0,1);
    else if constexpr(is_same_v<T,double>)
        ret=LAPACKE_dgeev(LAPACK_COL_MAJOR,'N','N',M.m,M.x,M.s,
            r.Get_Array_Pointer(),i.Get_Array_Pointer(),0,1,0,1);
    else PHYSBAM_FATAL_ERROR("LAPACK only defined for float/double types");

    PHYSBAM_ASSERT(!ret);
    for(int k=0;k<eigenvalues.m;k++)
        eigenvalues(k)=std::complex<T>(r(k),i(k));
#else
    PHYSBAM_FATAL_ERROR("This requies USE_LAPACK");
#endif
}

//#####################################################################
// CBLAS wrapper definition
//#####################################################################
template<class T>
void Matrix_Vector_Multiply_(MATRIX_VIEW<const T> M,ARRAY_VIEW<const T> A,
    ARRAY_VIEW<T> sol,T alpha,T beta,bool use_transpose)
{
#ifdef USE_LAPACK
    assert(M.Rows()==sol.m);
    assert(M.Columns()==A.m);
    CBLAS_TRANSPOSE transA=(use_transpose?CblasTrans:CblasNoTrans);

    if constexpr(is_same_v<T,float>)
        cblas_sgemv(CblasColMajor,transA,M.m,M.n,alpha,
            M.x, M.m,A.base_pointer,1,beta,sol.base_pointer,1);
    else if constexpr(is_same_v<T,double>)
        cblas_dgemv(CblasColMajor,transA,M.m,M.n,alpha,
            M.x, M.m,A.base_pointer,1,beta,sol.base_pointer,1);
    else PHYSBAM_FATAL_ERROR("LAPACK only defined for float/double types");

#else
    PHYSBAM_FATAL_ERROR("This requies USE_LAPACK");
#endif
}

template<class T>
void Symmetric_Matrix_Vector_Multiply_(MATRIX_VIEW<const T> M,
    ARRAY_VIEW<const T> A,ARRAY_VIEW<T> sol,T alpha,T beta)
{
#ifdef USE_LAPACK
    assert(M.Rows()==sol.m);
    assert(M.Columns()==A.m);
    assert(M.Columns()==M.Rows());

    if constexpr(is_same_v<T,float>)
        cblas_ssymv(CblasColMajor,CblasUpper,M.m,alpha,M.x, M.m,
            A.base_pointer,1,beta,sol.base_pointer,1);
    else if constexpr(is_same_v<T,double>)
        cblas_dsymv(CblasColMajor,CblasUpper,M.m,alpha,M.x, M.m,
            A.base_pointer,1,beta,sol.base_pointer,1);
    else PHYSBAM_FATAL_ERROR("LAPACK only defined for float/double types");

#else
    PHYSBAM_FATAL_ERROR("This requies USE_LAPACK");
#endif
}

//Destory A
template<class T>
void Triangular_Matrix_Vector_Multiply_(MATRIX_VIEW<const T> M,
    ARRAY_VIEW<T> A,bool use_upper,bool use_transpose,
    bool unit_diagonal)
{
#ifdef USE_LAPACK
    assert(M.Columns()==A.m);
    assert(M.Columns()==M.Rows());
    CBLAS_UPLO uplo=(use_upper?CblasUpper:CblasLower);
    CBLAS_TRANSPOSE transA=(use_transpose?CblasTrans:CblasNoTrans);
    CBLAS_DIAG diag=(unit_diagonal?CblasUnit:CblasNonUnit);

    if constexpr(is_same_v<T,float>)
        cblas_strmv(CblasColMajor,uplo,transA,diag,M.m,
            M.x, M.m,A.base_pointer,1);
    else if constexpr(is_same_v<T,double>)
        cblas_dtrmv(CblasColMajor,uplo,transA,diag,M.m,
            M.x, M.m,A.base_pointer,1);
    else PHYSBAM_FATAL_ERROR("LAPACK only defined for float/double types");

#else
    PHYSBAM_FATAL_ERROR("This requies USE_LAPACK");
#endif
}

//Destory A
template<class T>
void Triangular_Matrix_Problem_Solver_(MATRIX_VIEW<const T> M,
    ARRAY_VIEW<T> A,bool use_upper,bool use_transpose,
    bool unit_diagonal)
{
#ifdef USE_LAPACK
    assert(M.Columns()==A.m);    
    assert(M.Columns()==M.Rows());
    CBLAS_UPLO uplo=(use_upper?CblasUpper:CblasLower);
    CBLAS_TRANSPOSE transA=(use_transpose?CblasTrans:CblasNoTrans);
    CBLAS_DIAG diag=(unit_diagonal?CblasUnit:CblasNonUnit);

    if constexpr(is_same_v<T,float>)
        cblas_strsv(CblasColMajor,uplo,transA,diag,M.m,
            M.x, M.m,A.base_pointer,1);
    else if constexpr(is_same_v<T,double>)
        cblas_dtrsv(CblasColMajor,uplo,transA,diag,M.m,
            M.x, M.m,A.base_pointer,1);
    else PHYSBAM_FATAL_ERROR("LAPACK only defined for float/double types");

#else
    PHYSBAM_FATAL_ERROR("This requies USE_LAPACK");
#endif
}

//Destory M
template<class T>
void Rank_One_Operation_(MATRIX_VIEW<T> M,
    ARRAY_VIEW<const T> A,ARRAY_VIEW<const T> B, T alpha)
{
#ifdef USE_LAPACK
    assert(M.Columns()==A.Size());    
    assert(M.Rows()==B.Size());

    if constexpr(is_same_v<T,float>)
        cblas_sger(CblasColMajor,A.m,B.m,alpha,A.base_pointer,1,
            B.base_pointer,1,M.x,M.m);
    else if constexpr(is_same_v<T,double>)
        cblas_dger(CblasColMajor,A.m,B.m,alpha,A.base_pointer,1,
            B.base_pointer,1,M.x,M.m);
    else PHYSBAM_FATAL_ERROR("LAPACK only defined for float/double types");

#else
    PHYSBAM_FATAL_ERROR("This requies USE_LAPACK");
#endif
}

//Destory M
template<class T>
void Symmetric_Rank_One_Operation_(MATRIX_VIEW<T> M,ARRAY_VIEW<const T> A,
    T alpha)
{
#ifdef USE_LAPACK
    assert(M.Columns()==A.Size());    
    assert(M.Rows()==M.Columns());

    if constexpr(is_same_v<T,float>)
        cblas_ssyr(CblasColMajor,CblasLower,A.m,alpha,A.base_pointer,1,M.x,M.m);
    else if constexpr(is_same_v<T,double>)
        cblas_dsyr(CblasColMajor,CblasLower,A.m,alpha,A.base_pointer,1,M.x,M.m);
    else PHYSBAM_FATAL_ERROR("LAPACK only defined for float/double types");

    for(int i=0;i<M.m;i++) for(int j=0;j<i;j++) M(j,i)=M(i,j);
#else
    PHYSBAM_FATAL_ERROR("This requies USE_LAPACK");
#endif
}
    
//Destory M
template<class T>
void Symmetric_Rank_Two_Operation_(MATRIX_VIEW<T> M,ARRAY_VIEW<const T> A,
    ARRAY_VIEW<const T> B,T alpha)
{
#ifdef USE_LAPACK
    assert(M.Columns()==A.Size());    
    assert(M.Rows()==M.Columns());

    if constexpr(is_same_v<T,float>)
        cblas_ssyr2(CblasColMajor,CblasLower,A.m,alpha,A.base_pointer,1,
            B.base_pointer,1,M.x,M.m);
    else if constexpr(is_same_v<T,double>)
        cblas_dsyr2(CblasColMajor,CblasLower,A.m,alpha,A.base_pointer,1,
            B.base_pointer,1,M.x,M.m);
    else PHYSBAM_FATAL_ERROR("LAPACK only defined for float/double types");

    for(int i=0;i<M.m;i++) for(int j=0;j<i;j++) M(j,i)=M(i,j);

#else
    PHYSBAM_FATAL_ERROR("This requies USE_LAPACK");
#endif
}

//Destory C
template<class T>
void Matrix_Matrix_Multiply_(MATRIX_VIEW<const T> A,MATRIX_VIEW<const T> B,
    MATRIX_VIEW<T> C,T alpha,T beta,bool A_use_transpose,
    bool B_use_transpose)
{
#ifdef USE_LAPACK
    assert(C.Rows()==(A_use_transpose?A.Columns():A.Rows()));
    assert(C.Columns()==(B_use_transpose?B.Rows():B.Columns()));
    CBLAS_TRANSPOSE transA=(A_use_transpose?CblasTrans:CblasNoTrans);
    CBLAS_TRANSPOSE transB=(B_use_transpose?CblasTrans:CblasNoTrans);
    int k=(A_use_transpose?A.m:A.n);

    if constexpr(is_same_v<T,float>)
        cblas_sgemm(CblasColMajor,transA,transB,C.m,C.n,k,alpha,
            A.x,A.m,B.x,B.m,beta,C.x,C.m);
    else if constexpr(is_same_v<T,double>)
        cblas_dgemm(CblasColMajor,transA,transB,C.m,C.n,k,alpha,
            A.x,A.m,B.x,B.m,beta,C.x,C.m);
    else PHYSBAM_FATAL_ERROR("LAPACK only defined for float/double types");

#else
    PHYSBAM_FATAL_ERROR("This requies USE_LAPACK");
#endif
}

//Destory C
template<class T>
void Symmetric_Matrix_Matrix_Multiply_(MATRIX_VIEW<const T> A,MATRIX_VIEW<const T> B,
    MATRIX_VIEW<T> C,T alpha,T beta,bool A_on_left_side)
{
#ifdef USE_LAPACK
    assert(C.Rows()==(A_on_left_side?A.Rows():B.Rows()));
    assert(C.Columns()==(A_on_left_side?B.Columns():A.Columns()));
    CBLAS_SIDE side=(A_on_left_side?CblasLeft:CblasRight);

    if constexpr(is_same_v<T,float>)
        cblas_ssymm(CblasColMajor,side,CblasLower,C.m,C.n,alpha,A.x,A.m,
            B.x,B.m,beta,C.x,C.m);
    else if constexpr(is_same_v<T,double>)
        cblas_dsymm(CblasColMajor,side,CblasLower,C.m,C.n,alpha,A.x,A.m,
            B.x,B.m,beta,C.x,C.m);
    else PHYSBAM_FATAL_ERROR("LAPACK only defined for float/double types");

#else
    PHYSBAM_FATAL_ERROR("This requies USE_LAPACK");
#endif
}

//Destory C
template<class T>
void Symmetric_Rank_K_Operation_(MATRIX_VIEW<const T> A,MATRIX_VIEW<T> C,
    T alpha,T beta,bool use_transpose)
{
#ifdef USE_LAPACK
    assert(C.Rows()==(use_transpose?A.Columns():A.Rows()));
    assert(C.Rows()==C.Columns());
    CBLAS_TRANSPOSE trans=(use_transpose?CblasTrans:CblasNoTrans);
    int k=(use_transpose?A.Rows():A.Columns());

    if constexpr(is_same_v<T,float>)
        cblas_ssyrk(CblasColMajor,CblasLower,trans,C.m,k,alpha,A.x,A.m,
            beta,C.x,C.m);
    else if constexpr(is_same_v<T,double>)
        cblas_dsyrk(CblasColMajor,CblasLower,trans,C.m,k,alpha,A.x,A.m,
            beta,C.x,C.m);
    else PHYSBAM_FATAL_ERROR("LAPACK only defined for float/double types");

    for(int i=0;i<C.m;i++) for(int j=0;j<i;j++) C(j,i)=C(i,j);
    
#else
    PHYSBAM_FATAL_ERROR("This requies USE_LAPACK");
#endif
}

//Destory C
template<class T>
void Symmetric_Rank_2K_Operation_(MATRIX_VIEW<const T> A,MATRIX_VIEW<const T> B,
    MATRIX_VIEW<T> C,T alpha,T beta,bool use_transpose)
{
#ifdef USE_LAPACK
    assert(C.Rows()==(use_transpose?A.Columns():A.Rows()));
    assert(C.Rows()==C.Columns());
    assert(A.Rows()==B.Rows()); assert(A.Columns()==B.Columns());
    CBLAS_TRANSPOSE trans=(use_transpose?CblasTrans:CblasNoTrans);
    int k=(use_transpose?A.Rows():A.Columns());

    if constexpr(is_same_v<T,float>)
        cblas_ssyr2k(CblasColMajor,CblasLower,trans,C.m,k,alpha,A.x,A.m,
            B.x,B.m,beta,C.x,C.m);
    else if constexpr(is_same_v<T,double>)
        cblas_dsyr2k(CblasColMajor,CblasLower,trans,C.m,k,alpha,A.x,A.m,
            B.x,B.m,beta,C.x,C.m);
    else PHYSBAM_FATAL_ERROR("LAPACK only defined for float/double types");

    for(int i=0;i<C.m;i++) for(int j=0;j<i;j++) C(j,i)=C(i,j);

    
#else
    PHYSBAM_FATAL_ERROR("This requies USE_LAPACK");
#endif
}

//Destory B
template<class T>
void Triangular_Matrix_Matrix_Multiply_(MATRIX_VIEW<const T> A,
    MATRIX_VIEW<T> B,T alpha, bool A_on_left_side,bool use_upper,
    bool use_transpose,bool unit_diagonal)
{
#ifdef USE_LAPACK
    assert(A.Columns()==(A_on_left_side?B.Rows():B.Columns()));
    assert(A.Columns()==A.Rows());
    CBLAS_SIDE side=(A_on_left_side?CblasLeft:CblasRight);
    CBLAS_UPLO uplo=(use_upper?CblasUpper:CblasLower);
    CBLAS_TRANSPOSE transA=(use_transpose?CblasTrans:CblasNoTrans);
    CBLAS_DIAG diag=(unit_diagonal?CblasUnit:CblasNonUnit);

    if constexpr(is_same_v<T,float>)
        cblas_strmm(CblasColMajor,side,uplo,transA,diag,B.m,B.n,alpha,
            A.x,A.m,B.x,B.m);
    else if constexpr(is_same_v<T,double>)
        cblas_dtrmm(CblasColMajor,side,uplo,transA,diag,B.m,B.n,alpha,
            A.x,A.m,B.x,B.m);
    else PHYSBAM_FATAL_ERROR("LAPACK only defined for float/double types");

#else
    PHYSBAM_FATAL_ERROR("This requies USE_LAPACK");
#endif
}

//Destory B
template<class T>
void Triangular_Matrix_Solver_Multiple_(MATRIX_VIEW<const T> A,
    MATRIX_VIEW<T> B,T alpha,bool A_on_left_side,bool use_upper,
    bool use_transpose,bool unit_diagonal)
{
#ifdef USE_LAPACK
    assert(A.Columns()==(A_on_left_side?B.Rows():B.Columns()));
    assert(A.Columns()==A.Rows());
    CBLAS_SIDE side=(A_on_left_side?CblasLeft:CblasRight);
    CBLAS_UPLO uplo=(use_upper?CblasUpper:CblasLower);
    CBLAS_TRANSPOSE transA=(use_transpose?CblasTrans:CblasNoTrans);
    CBLAS_DIAG diag=(unit_diagonal?CblasUnit:CblasNonUnit);

    if constexpr(is_same_v<T,float>)
        cblas_strsm(CblasColMajor,side,uplo,transA,diag,B.m,B.n,alpha,
            A.x,A.m,B.x,B.m);
    else if constexpr(is_same_v<T,double>)
        cblas_dtrsm(CblasColMajor,side,uplo,transA,diag,B.m,B.n,alpha,
            A.x,A.m,B.x,B.m);
    else PHYSBAM_FATAL_ERROR("LAPACK only defined for float/double types");

#else
    PHYSBAM_FATAL_ERROR("This requies USE_LAPACK");
#endif
}

//#####################################################################

template void Compute_Eigenvalues_Raw_(MATRIX_VIEW<double>,
    ARRAY_VIEW<std::complex<double>>);
template void Compute_Eigenvalues_Raw_(MATRIX_VIEW<float>,
    ARRAY_VIEW<std::complex<float>>);
template void Least_Squares_Solve_Raw_(MATRIX_VIEW<double>,
    MATRIX_VIEW<double const>,MATRIX_VIEW<double>);
template void Least_Squares_Solve_Raw_(MATRIX_VIEW<float>,
    MATRIX_VIEW<float const>,MATRIX_VIEW<float>);
template void Matrix_Vector_Multiply_(MATRIX_VIEW<double const>,
    ARRAY_VIEW<double const>,ARRAY_VIEW<double>,double,double,bool);
template void Matrix_Vector_Multiply_(MATRIX_VIEW<float const>,
    ARRAY_VIEW<float const>,ARRAY_VIEW<float>,float,float,bool);
template void Symmetric_Matrix_Vector_Multiply_(MATRIX_VIEW<double const>,
    ARRAY_VIEW<double const>,ARRAY_VIEW<double>,double,double);
template void Symmetric_Matrix_Vector_Multiply_(MATRIX_VIEW<float const>,
    ARRAY_VIEW<float const>,ARRAY_VIEW<float>,float,float);
template void Triangular_Matrix_Vector_Multiply_(MATRIX_VIEW<double const>,
    ARRAY_VIEW<double>,bool,bool,bool);
template void Triangular_Matrix_Vector_Multiply_(MATRIX_VIEW<float const>,
    ARRAY_VIEW<float>,bool,bool,bool);
template void Triangular_Matrix_Problem_Solver_(MATRIX_VIEW<double const>,
    ARRAY_VIEW<double>,bool,bool,bool);
template void Triangular_Matrix_Problem_Solver_(MATRIX_VIEW<float const>,
    ARRAY_VIEW<float>,bool,bool,bool);
template void Rank_One_Operation_(MATRIX_VIEW<double>,
    ARRAY_VIEW<double const>,ARRAY_VIEW<double const>,double);
template void Rank_One_Operation_(MATRIX_VIEW<float>,
    ARRAY_VIEW<float const>,ARRAY_VIEW<float const>,float);
template void Symmetric_Rank_One_Operation_(MATRIX_VIEW<double>,
    ARRAY_VIEW<double const>,double);
template void Symmetric_Rank_One_Operation_(MATRIX_VIEW<float>,
    ARRAY_VIEW<float const>,float);
template void Symmetric_Rank_Two_Operation_(MATRIX_VIEW<double>,
    ARRAY_VIEW<double const>,ARRAY_VIEW<double const>,double);
template void Symmetric_Rank_Two_Operation_(MATRIX_VIEW<float>,
    ARRAY_VIEW<float const>,ARRAY_VIEW<float const>,float);
template void Matrix_Matrix_Multiply_(MATRIX_VIEW<double const>,
    MATRIX_VIEW<double const>,MATRIX_VIEW<double>,double,double,bool,bool);
template void Matrix_Matrix_Multiply_(MATRIX_VIEW<float const>,
    MATRIX_VIEW<float const>,MATRIX_VIEW<float>,float,float,bool,bool);
template void Symmetric_Matrix_Matrix_Multiply_(MATRIX_VIEW<double const>,
    MATRIX_VIEW<double const>,MATRIX_VIEW<double>,double,double,bool);
template void Symmetric_Matrix_Matrix_Multiply_(MATRIX_VIEW<float const>,
    MATRIX_VIEW<float const>,MATRIX_VIEW<float>,float,float,bool);
template void Symmetric_Rank_K_Operation_(MATRIX_VIEW<double const>,
    MATRIX_VIEW<double>,double,double,bool);
template void Symmetric_Rank_K_Operation_(MATRIX_VIEW<float const>,
    MATRIX_VIEW<float>,float,float,bool);
template void Symmetric_Rank_2K_Operation_(MATRIX_VIEW<double const>,
    MATRIX_VIEW<double const>,MATRIX_VIEW<double>,double,double,bool);
template void Symmetric_Rank_2K_Operation_(MATRIX_VIEW<float const>,
    MATRIX_VIEW<float const>,MATRIX_VIEW<float>,float,float,bool);
template void Triangular_Matrix_Matrix_Multiply_(MATRIX_VIEW<double const>,
    MATRIX_VIEW<double>,double,bool,bool,bool,bool);
template void Triangular_Matrix_Matrix_Multiply_(MATRIX_VIEW<float const>,
    MATRIX_VIEW<float>,float,bool,bool,bool,bool);
template void Triangular_Matrix_Solver_Multiple_(MATRIX_VIEW<double const>,
    MATRIX_VIEW<double>,double,bool,bool,bool,bool);
template void Triangular_Matrix_Solver_Multiple_(MATRIX_VIEW<float const>,
    MATRIX_VIEW<float>,float,bool,bool,bool,bool);
template void Least_Squares_Solve_Raw_(MATRIX_VIEW<std::complex<float> >,
    MATRIX_VIEW<std::complex<float> const>,
    MATRIX_VIEW<std::complex<float>>);
template void Least_Squares_Solve_Raw_(MATRIX_VIEW<std::complex<double> >,
    MATRIX_VIEW<std::complex<double> const>,
    MATRIX_VIEW<std::complex<double>>);
}
