Matrix.hpp

// Reaktoro is a unified framework for modeling chemically reactive systems.
//
// Copyright © 2014-2024 Allan Leal
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with this library. If not, see <http://www.gnu.org/licenses/>.
 
#pragma once
 
// Eigen includes
#include <Eigen/Core>
 
// autodiff includes
#include <autodiff/forward/real/eigen.hpp>
 
// Reaktoro includes
#include <Reaktoro/Common/Index.hpp>
#include <Reaktoro/Common/Real.hpp>
#include <Reaktoro/Common/TraitsUtils.hpp>
#include <Reaktoro/Common/TypeOp.hpp>
 
namespace Reaktoro {
 
template<typename T>
using VectorX = Eigen::Matrix<T, -1, 1, 0, -1, 1>;
 
template<typename T>
using MatrixX = Eigen::Matrix<T, -1, -1, 0, -1, -1>;
 
template<typename T>
using ArrayX = Eigen::Array<T, -1, 1, 0, -1, 1>;
 
template<typename T>
using ArrayXX = Eigen::Array<T, -1, -1, 0, -1, -1>;
 
template<typename T>
using RowVectorX = Eigen::Matrix<T, 1, -1, 0, 1, -1>;
 
//---------------------------------------------------------------------------------------------------------------------
// == VECTOR TYPE ALIASES ==
//---------------------------------------------------------------------------------------------------------------------
 
using VectorXr                = autodiff::VectorXreal;                               
using VectorXrRef             = Eigen::Ref<VectorXr>;                                
using VectorXrConstRef        = Eigen::Ref<const VectorXr>;                          
using VectorXrMap             = Eigen::Map<VectorXr>;                                
using VectorXrConstMap        = Eigen::Map<const VectorXr>;                          
using VectorXrStridedRef      = Eigen::Ref<VectorXr, 0, Eigen::InnerStride<>>;       
using VectorXrStridedConstRef = Eigen::Ref<const VectorXr, 0, Eigen::InnerStride<>>; 
 
using VectorXl                = VectorX<Eigen::Index>;                               
using VectorXlRef             = Eigen::Ref<VectorXl>;                                
using VectorXlConstRef        = Eigen::Ref<const VectorXl>;                          
using VectorXlMap             = Eigen::Map<VectorXl>;                                
using VectorXlConstMap        = Eigen::Map<const VectorXl>;                          
using VectorXlStridedRef      = Eigen::Ref<VectorXl, 0, Eigen::InnerStride<>>;       
using VectorXlStridedConstRef = Eigen::Ref<const VectorXl, 0, Eigen::InnerStride<>>; 
 
using VectorXd                = Eigen::VectorXd;                                     
using VectorXdRef             = Eigen::Ref<VectorXd>;                                
using VectorXdConstRef        = Eigen::Ref<const VectorXd>;                          
using VectorXdMap             = Eigen::Map<VectorXd>;                                
using VectorXdConstMap        = Eigen::Map<const VectorXd>;                          
using VectorXdStridedRef      = Eigen::Ref<VectorXd, 0, Eigen::InnerStride<>>;       
using VectorXdStridedConstRef = Eigen::Ref<const VectorXd, 0, Eigen::InnerStride<>>; 
 
 
//---------------------------------------------------------------------------------------------------------------------
// == ARRAY TYPE ALIASES ==
//---------------------------------------------------------------------------------------------------------------------
 
using ArrayXr                 = autodiff::ArrayXreal;                                
using ArrayXrRef              = Eigen::Ref<ArrayXr>;                                 
using ArrayXrConstRef         = Eigen::Ref<const ArrayXr>;                           
using ArrayXrMap              = Eigen::Map<ArrayXr>;                                 
using ArrayXrConstMap         = Eigen::Map<const ArrayXr>;                           
using ArrayXrStridedRef       = Eigen::Ref<ArrayXr, 0, Eigen::InnerStride<>>;        
using ArrayXrStridedConstRef  = Eigen::Ref<const ArrayXr, 0, Eigen::InnerStride<>>;  
 
using ArrayXl                 = ArrayX<Eigen::Index>;                                
using ArrayXlRef              = Eigen::Ref<ArrayXl>;                                 
using ArrayXlConstRef         = Eigen::Ref<const ArrayXl>;                           
using ArrayXlMap              = Eigen::Map<ArrayXl>;                                 
using ArrayXlConstMap         = Eigen::Map<const ArrayXl>;                           
using ArrayXlStridedRef       = Eigen::Ref<ArrayXl, 0, Eigen::InnerStride<>>;        
using ArrayXlStridedConstRef  = Eigen::Ref<const ArrayXl, 0, Eigen::InnerStride<>>;  
 
using ArrayXd                 = Eigen::ArrayXd;                                      
using ArrayXdRef              = Eigen::Ref<ArrayXd>;                                 
using ArrayXdConstRef         = Eigen::Ref<const ArrayXd>;                           
using ArrayXdMap              = Eigen::Map<ArrayXd>;                                 
using ArrayXdConstMap         = Eigen::Map<const ArrayXd>;                           
using ArrayXdStridedRef       = Eigen::Ref<ArrayXd, 0, Eigen::InnerStride<>>;        
using ArrayXdStridedConstRef  = Eigen::Ref<const ArrayXd, 0, Eigen::InnerStride<>>;  
 
using ArrayXXr                = autodiff::ArrayXXreal;                               
using ArrayXXrRef             = Eigen::Ref<ArrayXXr>;                                
using ArrayXXrConstRef        = Eigen::Ref<const ArrayXXr>;                          
using ArrayXXrMap             = Eigen::Map<ArrayXXr>;                                
using ArrayXXrConstMap        = Eigen::Map<const ArrayXXr>;                          
using ArrayXXrStridedRef      = Eigen::Ref<ArrayXXr, 0, Eigen::InnerStride<>>;       
using ArrayXXrStridedConstRef = Eigen::Ref<const ArrayXXr, 0, Eigen::InnerStride<>>; 
 
using ArrayXXd                = Eigen::ArrayXXd;                                     
using ArrayXXdRef             = Eigen::Ref<ArrayXXd>;                                
using ArrayXXdConstRef        = Eigen::Ref<const ArrayXXd>;                          
using ArrayXXdMap             = Eigen::Map<ArrayXXd>;                                
using ArrayXXdConstMap        = Eigen::Map<const ArrayXXd>;                          
using ArrayXXdStridedRef      = Eigen::Ref<ArrayXXd, 0, Eigen::InnerStride<>>;       
using ArrayXXdStridedConstRef = Eigen::Ref<const ArrayXXd, 0, Eigen::InnerStride<>>; 
 
//---------------------------------------------------------------------------------------------------------------------
// == MATRIX TYPE ALIASES ==
//---------------------------------------------------------------------------------------------------------------------
 
using MatrixXr                = autodiff::MatrixXreal;      
using MatrixXrRef             = Eigen::Ref<MatrixXr>;       
using MatrixXrConstRef        = Eigen::Ref<const MatrixXr>; 
using MatrixXrMap             = Eigen::Map<MatrixXr>;       
using MatrixXrConstMap        = Eigen::Map<const MatrixXr>; 
 
using MatrixXd                = Eigen::MatrixXd;            
using MatrixXdRef             = Eigen::Ref<MatrixXd>;       
using MatrixXdConstRef        = Eigen::Ref<const MatrixXd>; 
using MatrixXdMap             = Eigen::Map<MatrixXd>;       
using MatrixXdConstMap        = Eigen::Map<const MatrixXd>; 
 
//---------------------------------------------------------------------------------------------------------------------
// == ROW VECTOR TYPE ALIASES ==
//---------------------------------------------------------------------------------------------------------------------
 
using RowVectorXr                = autodiff::RowVectorXreal;      
using RowVectorXrRef             = Eigen::Ref<RowVectorXr>;       
using RowVectorXrConstRef        = Eigen::Ref<const RowVectorXr>; 
using RowVectorXrMap             = Eigen::Map<RowVectorXr>;       
using RowVectorXrConstMap        = Eigen::Map<const RowVectorXr>; 
 
using RowVectorXl                = RowVectorX<Eigen::Index>;      
using RowVectorXlRef             = Eigen::Ref<RowVectorXl>;       
using RowVectorXlConstRef        = Eigen::Ref<const RowVectorXl>; 
using RowVectorXlMap             = Eigen::Map<RowVectorXl>;       
using RowVectorXlConstMap        = Eigen::Map<const RowVectorXl>; 
 
using RowVectorXd                = Eigen::RowVectorXd;            
using RowVectorXdRef             = Eigen::Ref<RowVectorXd>;       
using RowVectorXdConstRef        = Eigen::Ref<const RowVectorXd>; 
using RowVectorXdMap             = Eigen::Map<RowVectorXd>;       
using RowVectorXdConstMap        = Eigen::Map<const RowVectorXd>; 
 
//---------------------------------------------------------------------------------------------------------------------
// == OTHER TYPE ALIASES ==
//---------------------------------------------------------------------------------------------------------------------
 
using PermutationMatrix = Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic>;
 
//---------------------------------------------------------------------------------------------------------------------
// == FUNCTION ALIASES ==
//---------------------------------------------------------------------------------------------------------------------
 
auto zeros(Index rows) -> decltype(VectorXd::Zero(rows));
 
auto ones(Index rows) -> decltype(VectorXd::Ones(rows));
 
auto random(Index rows) -> decltype(VectorXd::Random(rows));
 
auto linspace(double start, double stop, Index rows) -> decltype(VectorXd::LinSpaced(rows, start, stop));
 
auto unit(Index rows, Index i) -> decltype(VectorXd::Unit(rows, i));
 
auto zeros(Index rows, Index cols) -> decltype(MatrixXd::Zero(rows, cols));
 
auto ones(Index rows, Index cols) -> decltype(MatrixXd::Ones(rows, cols));
 
auto random(Index rows, Index cols) -> decltype(MatrixXd::Random(rows, cols));
 
auto identity(Index rows, Index cols) -> decltype(MatrixXd::Identity(rows, cols));
 
template<typename Derived>
auto rows(Eigen::MatrixBase<Derived>& mat, Index start, Index num) -> decltype(mat.middleRows(start, num));
 
template<typename Derived>
auto rows(const Eigen::MatrixBase<Derived>& mat, Index start, Index num) -> decltype(mat.middleRows(start, num));
 
template<typename Derived, typename Indices>
auto rows(Eigen::MatrixBase<Derived>& mat, const Indices& irows) -> decltype(mat(irows, Eigen::all));
 
template<typename Derived, typename Indices>
auto rows(const Eigen::MatrixBase<Derived>& mat, const Indices& irows) -> decltype(mat(irows, Eigen::all));
 
template<typename Derived>
auto cols(Eigen::MatrixBase<Derived>& mat, Index start, Index num) -> decltype(mat.middleCols(start, num));
 
template<typename Derived>
auto cols(const Eigen::MatrixBase<Derived>& mat, Index start, Index num) -> decltype(mat.middleCols(start, num));
 
template<typename Derived, typename Indices>
auto cols(Eigen::MatrixBase<Derived>& mat, const Indices& icols) -> decltype(mat(Eigen::all, icols));
 
template<typename Derived, typename Indices>
auto cols(const Eigen::MatrixBase<Derived>& mat, const Indices& icols) -> decltype(mat(Eigen::all, icols));
 
template<typename Derived>
auto segment(Eigen::MatrixBase<Derived>& vec, Index irow, Index nrows) -> decltype(vec.segment(irow, nrows));
 
template<typename Derived>
auto segment(const Eigen::MatrixBase<Derived>& vec, Index irow, Index nrows) -> decltype(vec.segment(irow, nrows));
 
template<typename Derived>
auto block(Eigen::MatrixBase<Derived>& mat, Index irow, Index icol, Index nrows, Index ncols) -> decltype(mat.block(irow, icol, nrows, ncols));
 
template<typename Derived>
auto block(const Eigen::MatrixBase<Derived>& mat, Index irow, Index icol, Index nrows, Index ncols) -> decltype(mat.block(irow, icol, nrows, ncols));
 
template<typename Derived, typename Indices>
auto submatrix(Eigen::MatrixBase<Derived>& mat, const Indices& irows, const Indices& icols) -> decltype(mat(irows, icols));
 
template<typename Derived, typename Indices>
auto submatrix(const Eigen::MatrixBase<Derived>& mat, const Indices& irows, const Indices& icols) -> decltype(mat(irows, icols));
 
template<typename Derived>
auto tr(Eigen::MatrixBase<Derived>& mat) -> decltype(mat.transpose());
 
template<typename Derived>
auto tr(const Eigen::MatrixBase<Derived>& mat) -> decltype(mat.transpose());
 
template<typename Derived>
auto inv(const Eigen::MatrixBase<Derived>& mat) -> decltype(mat.cwiseInverse());
 
template<typename Derived>
auto diag(const Eigen::MatrixBase<Derived>& vec) -> decltype(vec.asDiagonal());
 
template<typename Derived>
auto diagonal(Eigen::MatrixBase<Derived>& mat) -> decltype(mat.diagonal());
 
template<typename Derived>
auto diagonal(const Eigen::MatrixBase<Derived>& mat) -> decltype(mat.diagonal());
 
template<int p, typename Derived>
auto norm(const Eigen::MatrixBase<Derived>& mat) -> double;
 
template<typename Derived>
auto norm(const Eigen::MatrixBase<Derived>& mat) -> double;
 
template<typename Derived>
auto norminf(const Eigen::MatrixBase<Derived>& mat) -> double;
 
template<typename Derived>
auto sum(const Eigen::DenseBase<Derived>& mat) -> typename Derived::Scalar;
 
template<typename DerivedLHS, typename DerivedRHS>
auto dot(const Eigen::MatrixBase<DerivedLHS>& lhs, const Eigen::MatrixBase<DerivedRHS>& rhs) -> decltype(lhs.dot(rhs));
 
template<typename Derived>
auto min(const Eigen::MatrixBase<Derived>& mat) -> decltype(mat.minCoeff());
 
template<typename DerivedLHS, typename DerivedRHS>
auto min(const Eigen::MatrixBase<DerivedLHS>& lhs, const Eigen::MatrixBase<DerivedRHS>& rhs) -> decltype(lhs.cwiseMin(rhs));
 
template<typename Derived>
auto max(const Eigen::MatrixBase<Derived>& mat) -> decltype(mat.maxCoeff());
 
template<typename DerivedLHS, typename DerivedRHS>
auto max(const Eigen::MatrixBase<DerivedLHS>& lhs, const Eigen::MatrixBase<DerivedRHS>& rhs) -> decltype(lhs.cwiseMax(rhs));
 
template<typename DerivedLHS, typename DerivedRHS>
auto operator%(const Eigen::MatrixBase<DerivedLHS>& lhs, const Eigen::MatrixBase<DerivedRHS>& rhs) -> decltype(lhs.cwiseProduct(rhs));
 
template<typename DerivedLHS, typename DerivedRHS>
auto operator/(const Eigen::MatrixBase<DerivedLHS>& lhs, const Eigen::MatrixBase<DerivedRHS>& rhs) -> decltype(lhs.cwiseQuotient(rhs));
 
template<typename Derived>
auto operator/(const typename Derived::Scalar& scalar, const Eigen::MatrixBase<Derived>& mat) -> decltype(scalar*mat.cwiseInverse());
 
template<typename Derived>
auto operator+(const typename Derived::Scalar& scalar, const Eigen::MatrixBase<Derived>& mat) -> decltype((scalar + mat.array()).matrix());
 
template<typename Derived>
auto operator+(const Eigen::MatrixBase<Derived>& mat, const typename Derived::Scalar& scalar) -> decltype((scalar + mat.array()).matrix());
 
template<typename Derived>
auto operator-(const typename Derived::Scalar& scalar, const Eigen::MatrixBase<Derived>& mat) -> decltype((scalar - mat.array()).matrix());
 
template<typename Derived>
auto operator-(const Eigen::MatrixBase<Derived>& mat, const typename Derived::Scalar& scalar) -> decltype((mat.array() - scalar).matrix());
 
template<typename Derived>
auto abs(const Eigen::MatrixBase<Derived>& mat) -> decltype(mat.cwiseAbs());
 
template<typename Derived>
auto sqrt(const Eigen::MatrixBase<Derived>& mat) -> decltype(mat.cwiseSqrt());
 
template<typename Derived>
auto pow(const Eigen::MatrixBase<Derived>& mat, double power) -> decltype(mat.array().pow(power));
 
template<typename Derived>
auto exp(const Eigen::MatrixBase<Derived>& mat) -> decltype(mat.array().exp().matrix());
 
template<typename Derived>
auto log(const Eigen::MatrixBase<Derived>& mat) -> decltype(mat.array().log().matrix());
 
template<typename Derived>
auto log10(const Eigen::MatrixBase<Derived>& mat) -> decltype(mat.array().log10().matrix());
 
} // namespace Reaktoro
 
//=========================================================================
// CODE BELOW NEEDED FOR MEMOIZATION TECHNIQUE INVOLVING EIGEN TYPES
//=========================================================================
namespace Reaktoro {
 
template<typename T>
struct MemoizationTraits;
 
template<typename Scalar, int Rows, int Cols, int Options, int MaxRows, int MaxCols>
struct MemoizationTraits<Eigen::Array<Scalar, Rows, Cols, Options, MaxRows, MaxCols>>
{
    using Type = Eigen::Array<Scalar, Rows, Cols, Options, MaxRows, MaxCols>;
 
    using CacheType = Eigen::Array<Scalar, Rows, Cols, Options, MaxRows, MaxCols>;
 
    static auto equal(const CacheType& a, const Type& b)
    {
        return (a == b).all();
    }
};
 
template<typename EigenType>
struct MemoizationTraits<Eigen::Ref<EigenType>>
{
    using Type = Eigen::Ref<EigenType>;
 
    using CacheType = typename MemoizationTraits<Decay<EigenType>>::CacheType;
 
    static auto equal(const CacheType& a, const Type& b)
    {
        return MemoizationTraits<CacheType>::equal(a, b);
    }
 
    static auto assign(const CacheType& a, const Type& b)
    {
        return MemoizationTraits<CacheType>::assign(a, b);
    }
};
 
} // namespace Reaktoro
 
//=========================================================================
// CODE BELOW NEEDED FOR TypeOp TRAITS INVOLVING EIGEN TYPES
//=========================================================================
namespace Reaktoro {
namespace detail {
 
template<typename Scalar, int Rows, int Cols, int Options, int MaxRows, int MaxCols>
struct TypeOpRef<Eigen::Matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols>>
{
    using type = Eigen::Ref<Eigen::Matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols>>;
};
 
template<typename Scalar, int Rows, int Cols, int Options, int MaxRows, int MaxCols>
struct TypeOpRef<Eigen::Array<Scalar, Rows, Cols, Options, MaxRows, MaxCols>>
{
    using type = Eigen::Ref<Eigen::Array<Scalar, Rows, Cols, Options, MaxRows, MaxCols>>;
};
 
template<typename Scalar, int Rows, int Cols, int Options, int MaxRows, int MaxCols>
struct TypeOpConstRef<Eigen::Matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols>>
{
    using type = Eigen::Ref<const Eigen::Matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols>>;
};
 
template<typename Scalar, int Rows, int Cols, int Options, int MaxRows, int MaxCols>
struct TypeOpConstRef<Eigen::Array<Scalar, Rows, Cols, Options, MaxRows, MaxCols>>
{
    using type = Eigen::Ref<const Eigen::Array<Scalar, Rows, Cols, Options, MaxRows, MaxCols>>;
};
 
} // namespace detail
} // namespace Reaktoro
 
namespace Reaktoro {
 
inline auto zeros(Index rows) -> decltype(VectorXd::Zero(rows))
{
    return VectorXd::Zero(rows);
}
 
inline auto ones(Index rows) -> decltype(VectorXd::Ones(rows))
{
    return VectorXd::Ones(rows);
}
 
inline auto constants(Index rows, double val) -> decltype(VectorXd::Constant(rows, val))
{
    return VectorXd::Constant(rows, val);
}
 
inline auto random(Index rows) -> decltype(VectorXd::Random(rows))
{
    return VectorXd::Random(rows);
}
 
inline auto linspace(double start, double stop, Index rows) -> decltype(VectorXd::LinSpaced(rows, start, stop))
{
    return VectorXd::LinSpaced(rows, start, stop);
}
 
inline auto unit(Index rows, Index i) -> decltype(VectorXd::Unit(rows, i))
{
    return VectorXd::Unit(rows, i);
}
 
inline auto unitcol(Index rows, Index i) -> decltype(unit(rows, i))
{
    return unit(rows, i);
}
 
inline auto unitrow(Index cols, Index i) -> decltype(RowVectorXd::Unit(cols, i))
{
    return RowVectorXd::Unit(cols, i);
}
 
inline auto zeros(Index rows, Index cols) -> decltype(MatrixXd::Zero(rows, cols))
{
    return MatrixXd::Zero(rows, cols);
}
 
inline auto ones(Index rows, Index cols) -> decltype(MatrixXd::Ones(rows, cols))
{
    return MatrixXd::Ones(rows, cols);
}
 
inline auto random(Index rows, Index cols) -> decltype(MatrixXd::Random(rows, cols))
{
    return MatrixXd::Random(rows, cols);
}
 
inline auto identity(Index rows, Index cols) -> decltype(MatrixXd::Identity(rows, cols))
{
    return MatrixXd::Identity(rows, cols);
}
 
template<typename Derived>
auto rows(Eigen::MatrixBase<Derived>& mat, Index start, Index num) -> decltype(mat.middleRows(start, num))
{
    return mat.middleRows(start, num);
}
 
template<typename Derived>
auto rows(const Eigen::MatrixBase<Derived>& mat, Index start, Index num) -> decltype(mat.middleRows(start, num))
{
    return mat.middleRows(start, num);
}
 
template<typename Derived, typename Indices>
auto rows(Eigen::MatrixBase<Derived>& mat, const Indices& irows) -> decltype(mat(irows, Eigen::all))
{
    return mat(irows, Eigen::all);
}
 
template<typename Derived, typename Indices>
auto rows(const Eigen::MatrixBase<Derived>& mat, const Indices& irows) -> decltype(mat(irows, Eigen::all))
{
    return mat(irows, Eigen::all);
}
 
template<typename Derived>
auto cols(Eigen::MatrixBase<Derived>& mat, Index start, Index num) -> decltype(mat.middleCols(start, num))
{
    return mat.middleCols(start, num);
}
 
template<typename Derived>
auto cols(const Eigen::MatrixBase<Derived>& mat, Index start, Index num) -> decltype(mat.middleCols(start, num))
{
    return mat.middleCols(start, num);
}
 
template<typename Derived, typename Indices>
auto cols(Eigen::MatrixBase<Derived>& mat, const Indices& icols) -> decltype(mat(Eigen::all, icols))
{
    return mat(Eigen::all, icols);
}
 
template<typename Derived, typename Indices>
auto cols(const Eigen::MatrixBase<Derived>& mat, const Indices& icols) -> decltype(mat(Eigen::all, icols))
{
    return mat(Eigen::all, icols);
}
 
template<typename Derived>
auto segment(Eigen::MatrixBase<Derived>& vec, Index irow, Index nrows) -> decltype(vec.segment(irow, nrows))
{
    return vec.segment(irow, nrows);
}
 
template<typename Derived>
auto segment(const Eigen::MatrixBase<Derived>& vec, Index irow, Index nrows) -> decltype(vec.segment(irow, nrows))
{
    return vec.segment(irow, nrows);
}
 
template<typename Derived>
auto block(Eigen::MatrixBase<Derived>& mat, Index irow, Index icol, Index nrows, Index ncols) -> decltype(mat.block(irow, icol, nrows, ncols))
{
    return mat.block(irow, icol, nrows, ncols);
}
 
template<typename Derived>
auto block(const Eigen::MatrixBase<Derived>& mat, Index irow, Index icol, Index nrows, Index ncols) -> decltype(mat.block(irow, icol, nrows, ncols))
{
    return mat.block(irow, icol, nrows, ncols);
}
 
template<typename Derived, typename Indices>
auto submatrix(Eigen::MatrixBase<Derived>& mat, const Indices& irows, const Indices& icols) -> decltype(mat(irows, icols))
{
    return mat(irows, icols);
}
 
template<typename Derived, typename Indices>
auto submatrix(const Eigen::MatrixBase<Derived>& mat, const Indices& irows, const Indices& icols) -> decltype(mat(irows, icols))
{
    return mat(irows, icols);
}
 
template<typename Derived>
auto tr(Eigen::MatrixBase<Derived>& mat) -> decltype(mat.transpose())
{
    return mat.transpose();
}
 
template<typename Derived>
auto tr(const Eigen::MatrixBase<Derived>& mat) -> decltype(mat.transpose())
{
    return mat.transpose();
}
 
template<typename Derived>
auto inv(const Eigen::MatrixBase<Derived>& mat) -> decltype(mat.cwiseInverse())
{
    return mat.cwiseInverse();
}
 
template<typename Derived>
auto diag(const Eigen::MatrixBase<Derived>& vec) -> decltype(vec.asDiagonal())
{
    return vec.asDiagonal();
}
 
template<typename Derived>
auto diagonal(Eigen::MatrixBase<Derived>& mat) -> decltype(mat.diagonal())
{
    return mat.diagonal();
}
 
template<typename Derived>
auto diagonal(const Eigen::MatrixBase<Derived>& mat) -> decltype(mat.diagonal())
{
    return mat.diagonal();
}
 
template<int p, typename Derived>
auto norm(const Eigen::MatrixBase<Derived>& mat) -> double
{
    return mat.template lpNorm<p>();
}
 
template<typename Derived>
auto norm(const Eigen::MatrixBase<Derived>& mat) -> double
{
    return mat.norm();
}
 
template<typename Derived>
auto norminf(const Eigen::MatrixBase<Derived>& mat) -> double
{
    return mat.template lpNorm<Eigen::Infinity>();
}
 
template<typename Derived>
auto sum(const Eigen::DenseBase<Derived>& mat) -> typename Derived::Scalar
{
    return mat.sum();
}
 
template<typename DerivedLHS, typename DerivedRHS>
auto dot(const Eigen::MatrixBase<DerivedLHS>& lhs, const Eigen::MatrixBase<DerivedRHS>& rhs) -> decltype(lhs.dot(rhs))
{
    return lhs.dot(rhs);
}
 
template<typename Derived>
auto min(const Eigen::MatrixBase<Derived>& mat) -> decltype(mat.minCoeff())
{
    return mat.minCoeff();
}
 
template<typename DerivedLHS, typename DerivedRHS>
auto min(const Eigen::MatrixBase<DerivedLHS>& lhs, const Eigen::MatrixBase<DerivedRHS>& rhs) -> decltype(lhs.cwiseMin(rhs))
{
    return lhs.cwiseMin(rhs);
}
 
template<typename Derived>
auto max(const Eigen::MatrixBase<Derived>& mat) -> decltype(mat.maxCoeff())
{
    return mat.maxCoeff();
}
 
template<typename DerivedLHS, typename DerivedRHS>
auto max(const Eigen::MatrixBase<DerivedLHS>& lhs, const Eigen::MatrixBase<DerivedRHS>& rhs) -> decltype(lhs.cwiseMax(rhs))
{
    return lhs.cwiseMax(rhs);
}
 
template<typename DerivedLHS, typename DerivedRHS>
auto operator%(const Eigen::MatrixBase<DerivedLHS>& lhs, const Eigen::MatrixBase<DerivedRHS>& rhs) -> decltype(lhs.cwiseProduct(rhs))
{
    return lhs.cwiseProduct(rhs);
}
 
template<typename DerivedLHS, typename DerivedRHS>
auto operator/(const Eigen::MatrixBase<DerivedLHS>& lhs, const Eigen::MatrixBase<DerivedRHS>& rhs) -> decltype(lhs.cwiseQuotient(rhs))
{
    return lhs.cwiseQuotient(rhs);
}
 
template<typename Derived>
auto operator/(const typename Derived::Scalar& scalar, const Eigen::MatrixBase<Derived>& mat) -> decltype(scalar*mat.cwiseInverse())
{
    return scalar*mat.cwiseInverse();
}
 
template<typename Derived>
auto operator+(const typename Derived::Scalar& scalar, const Eigen::MatrixBase<Derived>& mat) -> decltype((scalar + mat.array()).matrix())
{
    return (scalar + mat.array()).matrix();
}
 
template<typename Derived>
auto operator+(const Eigen::MatrixBase<Derived>& mat, const typename Derived::Scalar& scalar) -> decltype((scalar + mat.array()).matrix())
{
    return (scalar + mat.array()).matrix();
}
 
template<typename Derived>
auto operator-(const typename Derived::Scalar& scalar, const Eigen::MatrixBase<Derived>& mat) -> decltype((scalar - mat.array()).matrix())
{
    return (scalar - mat.array()).matrix();
}
 
template<typename Derived>
auto operator-(const Eigen::MatrixBase<Derived>& mat, const typename Derived::Scalar& scalar) -> decltype((mat.array() - scalar).matrix())
{
    return (mat.array() - scalar).matrix();
}
 
template<typename Derived>
auto abs(const Eigen::MatrixBase<Derived>& mat) -> decltype(mat.cwiseAbs())
{
    return mat.cwiseAbs();
}
 
template<typename Derived>
auto sqrt(const Eigen::MatrixBase<Derived>& mat) -> decltype(mat.cwiseSqrt())
{
    return mat.cwiseSqrt();
}
 
template<typename Derived>
auto pow(const Eigen::MatrixBase<Derived>& mat, double power) -> decltype(mat.array().pow(power).matrix())
{
    return mat.array().pow(power).matrix();
}
 
template<typename Derived>
auto exp(const Eigen::MatrixBase<Derived>& mat) -> decltype(mat.array().exp().matrix())
{
    return mat.array().exp().matrix();
}
 
template<typename Derived>
auto log(const Eigen::MatrixBase<Derived>& mat) -> decltype(mat.array().log().matrix())
{
    return mat.array().log().matrix();
}
 
template<typename Derived>
auto log10(const Eigen::MatrixBase<Derived>& mat) -> decltype(mat.array().log10().matrix())
{
    return mat.array().log10().matrix();
}
 
} // namespace Reaktoro