//#####################################################################
// Copyright 2019, Craig Schroeder.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class LAPACK
//#####################################################################
#ifndef __LAPACK__
#define __LAPACK__
#include <Core/Arrays/ARRAY.h>
#include <Core/Arrays/ARRAY_VIEW.h>
#include <Core/Matrices/MATRIX_BASE.h>
#include <Core/Matrices/MATRIX_VIEW.h>
#include <complex>
#ifdef USE_LAPACK
#define lapack_complex_double std::complex<double>
#define lapack_complex_float std::complex<float>
extern "C"{
#include <cblas.h>
#include <lapacke.h>
}
#endif
namespace PhysBAM{

//#####################################################################
// Forward Linear_Solve
//#####################################################################

// Linear_Solve_Raw and Linear_Solve require square and non-singular M.

// Destroys M
template<class T>
void Linear_Solve_Raw_(MATRIX_VIEW<T> M,MATRIX_VIEW<const T> rhs,MATRIX_VIEW<T> sol);

template<class T,class T_MATRIX1,class T_MATRIX2,class T_MATRIX3>
void Linear_Solve_Raw(MATRIX_BASE<T,T_MATRIX1>& M,
    const MATRIX_BASE<T,T_MATRIX2>& rhs,MATRIX_BASE<T,T_MATRIX3>& sol)
{
    Linear_Solve_Raw_<T>(M.Derived(),rhs.Derived(),sol.Derived());
}

template<class T,class T_MATRIX1,class T_MATRIX2,class T_MATRIX3>
void Linear_Solve(const MATRIX_BASE<T,T_MATRIX1>& M,
    const MATRIX_BASE<T,T_MATRIX2>& rhs,MATRIX_BASE<T,T_MATRIX3>& sol)
{
    MAT_TYPE<T,mat_gen,T_MATRIX1::fixed_m,T_MATRIX1::fixed_n> MM(M.Derived());
    Linear_Solve_Raw_<T>(MM,rhs,sol);
}

// Destroys M
template<class T,class T_MATRIX,class T_ARRAY1,class T_ARRAY2>
void Linear_Solve_Raw(MATRIX_BASE<T,T_MATRIX>& M,
    const ARRAY_BASE<T,T_ARRAY1>& rhs,ARRAY_BASE<T,T_ARRAY2>& sol)
{
    MATRIX_VIEW<const T> r(ARRAY_VIEW<T>(rhs.Derived()));
    MATRIX_VIEW<T> s(ARRAY_VIEW<T>(sol.Derived()));
    Linear_Solve_Raw_<T>(M.Derived(),r,s);
}

// Preserves M
template<class T,class T_MATRIX,class T_ARRAY1,class T_ARRAY2>
void Linear_Solve(const MATRIX_BASE<T,T_MATRIX>& M,
    const ARRAY_BASE<T,T_ARRAY1>& rhs,ARRAY_BASE<T,T_ARRAY2>& sol)
{
    MAT_TYPE<T,mat_gen,T_MATRIX::fixed_m,T_MATRIX::fixed_n> MM(M.Derived());
    MATRIX_VIEW<const T> r(ARRAY_VIEW<T>(rhs.Derived()));
    MATRIX_VIEW<T> s(ARRAY_VIEW<T>(sol.Derived()));
    Linear_Solve_Raw_<T>(MM,r,s);
}

//#####################################################################
// Forward Least_Squares_Solve
//#####################################################################

// Least_Squares_Solve_Raw and Least_Squares_Solve work with non-square and singular M.

// Destroys M
template<class T>
void Least_Squares_Solve_Raw_(MATRIX_VIEW<T> M,MATRIX_VIEW<const T> rhs,MATRIX_VIEW<T> sol);

template<class T,class T_MATRIX1,class T_MATRIX2,class T_MATRIX3>
void Least_Squares_Solve_Raw(MATRIX_BASE<T,T_MATRIX1>& M,
    const MATRIX_BASE<T,T_MATRIX2>& rhs,MATRIX_BASE<T,T_MATRIX3>& sol)
{
    Least_Squares_Solve_Raw_<T>(M.Derived(),rhs.Derived(),sol.Derived());
}

template<class T,class T_MATRIX1,class T_MATRIX2,class T_MATRIX3>
void Least_Squares_Solve(const MATRIX_BASE<T,T_MATRIX1>& M,
    const MATRIX_BASE<T,T_MATRIX2>& rhs,MATRIX_BASE<T,T_MATRIX3>& sol)
{
    MAT_TYPE<T,mat_gen,T_MATRIX1::fixed_m,T_MATRIX1::fixed_n> MM(M.Derived());
    Least_Squares_Solve_Raw_<T>(MM,rhs,sol);
}

// Destroys M
template<class T,class T_MATRIX,class T_ARRAY1,class T_ARRAY2>
void Least_Squares_Solve_Raw(MATRIX_BASE<T,T_MATRIX>& M,
    const ARRAY_BASE<T,T_ARRAY1>& rhs,ARRAY_BASE<T,T_ARRAY2>& sol)
{
    MATRIX_VIEW<const T> r(ARRAY_VIEW<const T>(rhs.Derived()));
    MATRIX_VIEW<T> s(ARRAY_VIEW<T>(sol.Derived()));
    Least_Squares_Solve_Raw_<T>(M.Derived(),r,s);
}

// Preserves M
template<class T,class T_MATRIX,class T_ARRAY1,class T_ARRAY2>
void Least_Squares_Solve(const MATRIX_BASE<T,T_MATRIX>& M,
    const ARRAY_BASE<T,T_ARRAY1>& rhs,ARRAY_BASE<T,T_ARRAY2>& sol)
{
    MAT_TYPE<T,mat_gen,T_MATRIX::fixed_m,T_MATRIX::fixed_n> MM(M.Derived());
    MATRIX_VIEW<const T> r(ARRAY_VIEW<const T>(rhs.Derived()));
    MATRIX_VIEW<T> s(ARRAY_VIEW<T>(sol.Derived()));
    Least_Squares_Solve_Raw_<T>(MM,r,s);
}

// Compute (complex) eigenvalues for non-symmetric real-valued or complex-valued
// matrix.  Destroys M
template<class T>
void Compute_Eigenvalues_Raw_(MATRIX_VIEW<T> M,ARRAY_VIEW<std::complex<T>> ev);

template<class T,class T_MATRIX,class T_ARRAY>
void Compute_Eigenvalues_Raw(MATRIX_BASE<T,T_MATRIX>& M,
    ARRAY_BASE<std::complex<T>,T_ARRAY>& ev)
{
    Compute_Eigenvalues_Raw_<T>(M.Derived(),ev.Derived());
}

template<class T,class T_MATRIX,class T_ARRAY>
void Compute_Eigenvalues(const MATRIX_BASE<T,T_MATRIX>& M,
    ARRAY_BASE<std::complex<T>,T_ARRAY>& ev)
{
    MAT_TYPE<T,mat_gen,T_MATRIX::fixed_m,T_MATRIX::fixed_n> MM(M.Derived());
    Compute_Eigenvalues_Raw_<T>(MM,ev.Derived());
}

//#####################################################################
// CBLAS forward declaration
//#####################################################################

// if use_transpose: sol=alpha*M^T*A+beta*sol
// else: sol=alpha*M*A+beta*sol
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);

template<class T,class T_MATRIX,class T_ARRAY1,class T_ARRAY2>
void Matrix_Vector_Multiply(const MATRIX_BASE<T,T_MATRIX>& M,
    const ARRAY_BASE<T,T_ARRAY1>& A,ARRAY_BASE<T,T_ARRAY2>& sol,
    T alpha,T beta,bool use_transpose)
{
    Matrix_Vector_Multiply_<T>(M.Derived(),A.Derived(),sol.Derived(),alpha,beta,use_transpose);
}

// Assume M is a symmetric matrix
// Computes sol=alpha*M*A+beta*sol
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);

template<class T,class T_MATRIX,class T_ARRAY1,class T_ARRAY2>
void Symmetric_Matrix_Vector_Multiply(const MATRIX_BASE<T,T_MATRIX>& M,
    const ARRAY_BASE<T,T_ARRAY1>& A,ARRAY_BASE<T,T_ARRAY2>& sol,
    T alpha,T beta)
{
    Symmetric_Matrix_Vector_Multiply_<T>(M.Derived(),A.Derived(),sol.Derived(),alpha,beta);
}

// if use_upper, only look at the upper triangular port of M.
// Otherwise, only look at the lower triangular port of M.
// if use_transpose, computes A=M*A.  Otherwise, computes A=M^T*A.
// if unit_diagonal, A is a unit triangular matrix.
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);

template<class T,class T_MATRIX,class T_ARRAY>
void Triangular_Matrix_Vector_Multiply(const MATRIX_BASE<T,T_MATRIX>& M,
    ARRAY_BASE<T,T_ARRAY>& A,bool use_upper,bool use_transpose,
    bool unit_diagonal)
{
    Triangular_Matrix_Vector_Multiply_<T>(M.Derived(),A.Derived(),use_upper,
        use_transpose,unit_diagonal);
}

// if use_upper, only look at the upper triangular port of M.
// Otherwise, only look at the lower triangular port of M.
// if use_transpose, solves A*x=b.  Otherwise, solves A^T*x=b.
// if unit_diagonal, A is a unit triangular matrix.
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);

template<class T,class T_MATRIX,class T_ARRAY>
void Triangular_Matrix_Problem_Solver(const MATRIX_BASE<T,T_MATRIX>& M,
    ARRAY_BASE<T,T_ARRAY>& A,bool use_upper,bool use_transpose,
    bool unit_diagonal)
{
    Triangular_Matrix_Problem_Solver_<T>(M.Derived(),A.Derived(),use_upper,
        use_transpose,unit_diagonal);
}

// Computes M+=alpha*A*B^T
template<class T>
void Rank_One_Operation_(MATRIX_VIEW<T> M,
    ARRAY_VIEW<const T> A,ARRAY_VIEW<const T> B,T alpha);

template<class T,class T_MATRIX,class T_ARRAY1,class T_ARRAY2>
void Rank_One_Operation(MATRIX_BASE<T,T_MATRIX>& M,
    const ARRAY_BASE<T,T_ARRAY1>& A,const ARRAY_BASE<T,T_ARRAY2>& B,
    T alpha)
{
    Rank_One_Operation_<T>(M.Derived(),A.Derived(),B.Derived(),alpha);
}

// Assume M is a symmetric matrix. Computes M+=alpha*A*A^T
template<class T>
void Symmetric_Rank_One_Operation_(MATRIX_VIEW<T> M,
    ARRAY_VIEW<const T> A,T alpha);

template<class T,class T_MATRIX,class T_ARRAY>
void Symmetric_Rank_One_Operation(MATRIX_BASE<T,T_MATRIX>& M,
    const ARRAY_BASE<T,T_ARRAY>& A,T alpha)
{
    Symmetric_Rank_One_Operation_<T>(M.Derived(),A.Derived(),alpha);
}

// Assume M is a symmetric matrix.
// Computes M=alpha*A*B^T+alpha*B*A^T+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);

template<class T,class T_MATRIX,class T_ARRAY1,class T_ARRAY2>
void Symmetric_Rank_Two_Operation(MATRIX_BASE<T,T_MATRIX>& M,
    const ARRAY_BASE<T,T_ARRAY1>& A,const ARRAY_BASE<T,T_ARRAY2>& B,
    T alpha)
{
    Symmetric_Rank_Two_Operation_<T>(M.Derived(),A.Derived(),B.Derived(),alpha);
}

// if A_use_transpose and B_use_transpose: C=alpha*A^T*B^T+beta*C
// if A_use_transpose and not B_use_transpose: C=alpha*A^T*B+beta*C
// if not A_use_transpose and B_use_transpose: C=alpha*A*B^T+beta*C
// if not A_use_transpose and not B_use_transpose: C=alpha*A*B+beta*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);

template<class T,class T_MATRIX1,class T_MATRIX2,class T_MATRIX3>
void Matrix_Matrix_Multiply(const MATRIX_BASE<T,T_MATRIX1>& A,
    const MATRIX_BASE<T,T_MATRIX2>& B,MATRIX_BASE<T,T_MATRIX3>& C,
    T alpha,T beta,bool A_use_transpose,bool B_use_transpose)
{
    Matrix_Matrix_Multiply_<T>(A.Derived(),B.Derived(),C.Derived(),alpha,beta,
        A_use_transpose,B_use_transpose);
}

// Assume A is a symmetric matrix, B and C are m by n matrices.
// If A_on_left_side, computes C=alpha*A*B+beta*C
// Else, computes C=alpha*B*A+beta*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);

template<class T,class T_MATRIX1,class T_MATRIX2,class T_MATRIX3>
void Symmetric_Matrix_Matrix_Multiply(const MATRIX_BASE<T,T_MATRIX1>& A,
    const MATRIX_BASE<T,T_MATRIX2>& B,MATRIX_BASE<T,T_MATRIX3>& C,
    T alpha,T beta,bool A_on_left_side)
{
    Symmetric_Matrix_Matrix_Multiply_<T>(A.Derived(),B.Derived(),C.Derived(),
        alpha,beta,A_on_left_side);
}

// Assume C is an n by n symmetric matrix.
// If use_transpose, A is an n by k matrix,
// computes C=alpha*A^T*A+beta*C.
// Else, A is a k by n matrix, computes C=alpha*A*A^T+beta*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);

template<class T,class T_MATRIX1,class T_MATRIX2>
void Symmetric_Rank_K_Operation(const MATRIX_BASE<T,T_MATRIX1>& A,
    MATRIX_BASE<T,T_MATRIX2>& C,T alpha,T beta,bool use_transpose)
{
    Symmetric_Rank_K_Operation_<T>(A.Derived(),C.Derived(),alpha,beta,use_transpose);
}

// Assume C is an n by n symmetric matrix.
// If use_transpose, A and B are k by n matrices,
// computes C=alpha*A^T*B+alpha*B^T*A+beta*C
// Else, A and B are n by k matrices,
// computes C=alpha*A*B^T+alpha*B*A^T+beta*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);

template<class T,class T_MATRIX1,class T_MATRIX2,class T_MATRIX3>
void Symmetric_Rank_2K_Operation(const MATRIX_BASE<T,T_MATRIX1>& A,
    const MATRIX_BASE<T,T_MATRIX2>& B,MATRIX_BASE<T,T_MATRIX3>& C,
    T alpha,T beta,bool use_transpose)
{
    Symmetric_Rank_2K_Operation_<T>(A.Derived(),B.Derived(),C.Derived(),alpha,beta,use_transpose);
}

// if A_use_transpose and A_on_left_side: B=alpha*A^T*B.
// if not A_use_transpose and A_on_left_side: B=alpha*A*B.
// if A_use_transpose and not A_on_left_side: B=alpha*B*A^T.
// if not A_use_transpose and not A_on_left_side: B=alpha*B*A.
// if use_upper, assume A is an upper triangular matrix,
// else, assume A is an lower triangular matrix.
// if unit_diagonal, assume A is an unit triangular matrix.
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 A_use_transpose,bool unit_diagonal);

template<class T,class T_MATRIX1,class T_MATRIX2>
void Triangular_Matrix_Matrix_Multiply(const MATRIX_BASE<T,T_MATRIX1>& A,
    MATRIX_BASE<T,T_MATRIX2>& B,T alpha,bool A_on_left_side,
    bool use_upper,bool use_transpose,bool unit_diagonal)
{
    Triangular_Matrix_Matrix_Multiply_<T>(A.Derived(),B.Derived(),alpha,A_on_left_side,
        use_upper,use_transpose,unit_diagonal);
}

// Assume A is an nxn triangular matrix, X and B are m by n matrices.
// if A_use_transpose and A_on_left_side: A^T*X=alpha*B.
// if not A_use_transpose and A_on_left_side: A*X=alpha*B.
// if A_use_transpose and not A_on_left_side: X*A^T=alpha*B.
// if not A_use_transpose and not A_on_left_side: X*A^T=alpha*B.
// if use_upper, assume A is an upper triangular matrix,
// else, assume A is an lower triangular matrix.
// if unit_diagonal, assume A is an unit triangular.
// B is overwritten by X on exit.
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);

template<class T,class T_MATRIX1,class T_MATRIX2>
void Triangular_Matrix_Solver_Multiple(const MATRIX_BASE<T,T_MATRIX1>& A,
    MATRIX_BASE<T,T_MATRIX2>& B,T alpha,bool A_on_left_side,
    bool use_upper,bool use_transpose,bool unit_diagonal)
{
    Triangular_Matrix_Solver_Multiple_<T>(A.Derived(),B.Derived(),alpha,A_on_left_side,
        use_upper,use_transpose,unit_diagonal);
}
}
#endif
