Solve.h

Go to the documentation of this file.

#pragma once

 
#include <Eigen/Core>
 
#include <EmrysFilter/Core/Multiply.h>
 
namespace Emrys::internal
{

inline constexpr int unrolledSolveSize = 8;

template <typename Result, typename Decomposition, typename Rhs>
[[nodiscard]] Result
transposedSolve(Decomposition const & decomposition, Eigen::MatrixBase<Rhs> const & rhs)
{
    constexpr int size = Decomposition::MatrixType::RowsAtCompileTime;
    if constexpr (size != Eigen::Dynamic && size <= unrolledSolveSize)
    {
        using Column = Eigen::Matrix<typename Rhs::Scalar, size, 1>;
        auto result = Result{};
        result.resize(rhs.cols(), rhs.rows());
        for (Eigen::Index j = 0; j < rhs.cols(); ++j)
        {
            Column const column = rhs.col(j);
            result.row(j) = decomposition.solve(column).transpose();
        }
        return result;
    }
    else
    {
        Result result = decomposition.solve(rhs).transpose();
        return result;
    }
}

template <typename Solver, typename Derived>
void
solveInPlace(Solver const & solver, Eigen::MatrixBase<Derived> & matrix)
{
    constexpr int size = Derived::RowsAtCompileTime;
    if constexpr (size != Eigen::Dynamic && size <= unrolledSolveSize && Derived::ColsAtCompileTime != 1)
    {
        using Column = Eigen::Matrix<typename Derived::Scalar, size, 1>;
        for (Eigen::Index j = 0; j < matrix.cols(); ++j)
        {
            Column column = matrix.col(j);
            solver.solveInPlace(column);
            matrix.col(j) = column;
        }
    }
    else
    {
        solver.solveInPlace(matrix);
    }
}

template <typename Result, typename Decomposition>
[[nodiscard]] Result
inverseFrom(Decomposition const & decomposition)
{
    using Matrix = typename Decomposition::MatrixType;
    constexpr int size = Matrix::RowsAtCompileTime;
    Eigen::Index const n = decomposition.rows();
    if constexpr (size != Eigen::Dynamic && size <= unrolledSolveSize)
    {
        using Column = Eigen::Matrix<typename Matrix::Scalar, size, 1>;
        auto result = Result{};
        result.resize(n, n);
        for (Eigen::Index j = 0; j < n; ++j)
        {
            result.col(j) = decomposition.solve(Column::Unit(n, j));
        }
        return result;
    }
    else
    {
        Result result = decomposition.solve(Result::Identity(n, n));
        return result;
    }
}
 
}  // namespace Emrys::internal