gsMatrixOp.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsLinearAlgebra.h>
#include <gsSolver/gsLinearOperator.h>

namespace gismo
{

// left here for debugging purposes
// template<typename T> struct is_ref { static const bool value = false; };
// template<typename T> struct is_ref<T&> { static const bool value = true; };

template <class MatrixType>
class gsMatrixOp GISMO_FINAL : public gsLinearOperator<typename MatrixType::Scalar>
{
    typedef memory::shared_ptr<MatrixType> MatrixPtr;
    typedef typename MatrixType::Nested NestedMatrix;

public:
    typedef typename MatrixType::Scalar T;

    typedef memory::shared_ptr<gsMatrixOp> Ptr;

    typedef memory::unique_ptr<gsMatrixOp> uPtr;

    gsMatrixOp(const MatrixType& mat)
    : m_mat(), m_expr(mat.derived())
    {
        //gsDebug<<typeid(m_expr).name()<<" Ref: "<<is_ref<NestedMatrix>::value<<"\n";
    }

    gsMatrixOp(MatrixPtr mat)
    : m_mat(give(mat)), m_expr(m_mat->derived())
    { }

    static uPtr make(const MatrixType& mat)
    { return uPtr( new gsMatrixOp(mat) ); }

    static uPtr make(MatrixPtr mat)
    { return uPtr( new gsMatrixOp(give(mat)) ); }

    void apply(const gsMatrix<T> & input, gsMatrix<T> & x) const
    { x.noalias() = m_expr * input; }

    index_t rows() const
    { return m_expr.rows(); }

    index_t cols() const
    { return m_expr.cols(); }

    NestedMatrix matrix() const
    { return m_expr; }

    MatrixPtr    matrixPtr() const {
        GISMO_ENSURE( m_mat, "A shared pointer is only available if it was provided to gsMatrixOp." );
        return m_mat;
    }

private:
    const MatrixPtr m_mat; 
    NestedMatrix   m_expr; 
};

template <class Derived>
typename gsMatrixOp<Derived>::uPtr makeMatrixOp(const gsEigen::EigenBase<Derived>& mat)
{
    return gsMatrixOp<Derived>::make(mat.derived());
}

template <class Derived>
typename gsMatrixOp<Derived>::uPtr makeMatrixOp(memory::shared_ptr<Derived> mat)
{
    return memory::make_unique(new gsMatrixOp<Derived>(give(mat)));
}

// We need an additional guide for the compiler to be able to work well with unique ptrs
template <class Derived>
typename gsMatrixOp<Derived>::uPtr makeMatrixOp(memory::unique_ptr<Derived> mat)
{
    return memory::make_unique(new gsMatrixOp<Derived>(memory::shared_ptr<Derived>(mat.release())));
}

template <class SolverType>
class gsSolverOp GISMO_FINAL : public gsLinearOperator<typename SolverType::Scalar>
{
public:
    typedef typename SolverType::Scalar T;

    typedef typename SolverType::MatrixType MatrixType;

    typedef memory::shared_ptr<gsSolverOp> Ptr;

    typedef memory::unique_ptr<gsSolverOp> uPtr;


    gsSolverOp(const MatrixType& mat)
    {
        GISMO_ASSERT(mat.rows() == mat.cols(), "Need square matrix");
        m_size = mat.rows();

        m_solver.compute(mat);
    }

    gsSolverOp(const memory::shared_ptr<MatrixType>& mat)
    {
        GISMO_ASSERT(mat->rows() == mat->cols(), "Need square matrix");
        m_size = mat->rows();

        m_solver.compute(*mat);
    }

    static uPtr make(const MatrixType& mat) { return memory::make_unique( new gsSolverOp(mat) ); }

    void apply(const gsMatrix<T> & input, gsMatrix<T> & x) const
    {
        x = m_solver.solve(input);
    }

    index_t rows() const { return m_size; }

    index_t cols() const { return m_size; }

    SolverType& solver()                { return m_solver; }

    const SolverType& solver() const    { return m_solver; }

private:
    SolverType m_solver;
    index_t m_size;
};


template <class T, int _Rows, int _Cols, int _Opt>
typename gsSolverOp< gsEigen::PartialPivLU< gsEigen::Matrix<T, _Rows, _Cols, _Opt> > >::uPtr  makePartialPivLUSolver(const gsMatrix<T, _Rows, _Cols, _Opt> & mat)
{
    return memory::make_unique( new gsSolverOp< gsEigen::PartialPivLU< gsEigen::Matrix<T, _Rows, _Cols, _Opt> > >(mat) );
}

template <class T, int _Rows, int _Cols, int _Opt>
typename gsSolverOp< gsEigen::PartialPivLU< gsEigen::Matrix<T, _Rows, _Cols, _Opt> > >::uPtr  makePartialPivLUSolver(const memory::shared_ptr< gsMatrix<T, _Rows, _Cols, _Opt> > & mat)
{
    return memory::make_unique( new gsSolverOp< gsEigen::PartialPivLU< gsEigen::Matrix<T, _Rows, _Cols, _Opt> > >(mat) );
}


template <class T, int _Rows, int _Cols, int _Opt>
typename gsSolverOp< gsEigen::FullPivLU< gsEigen::Matrix<T, _Rows, _Cols, _Opt> > >::uPtr  makeFullPivLUSolver(const gsMatrix<T, _Rows, _Cols, _Opt> & mat)
{
    return memory::make_unique( new gsSolverOp< gsEigen::FullPivLU< gsEigen::Matrix<T, _Rows, _Cols, _Opt> > >(mat) );
}

template <class T, int _Rows, int _Cols, int _Opt>
typename gsSolverOp< gsEigen::FullPivLU< gsEigen::Matrix<T, _Rows, _Cols, _Opt> > >::uPtr  makeFullPivLUSolver(const memory::shared_ptr< gsMatrix<T, _Rows, _Cols, _Opt> > & mat)
{
    return memory::make_unique( new gsSolverOp< gsEigen::FullPivLU< gsEigen::Matrix<T, _Rows, _Cols, _Opt> > >(mat) );
}


template <class T, int _Rows, int _Cols, int _Opt>
typename gsSolverOp< gsEigen::LDLT< gsEigen::Matrix<T, _Rows, _Cols, _Opt> > >::uPtr  makeCholeskySolver(const gsMatrix<T, _Rows, _Cols, _Opt> & mat)
{
    return memory::make_unique( new gsSolverOp< gsEigen::LDLT< gsEigen::Matrix<T, _Rows, _Cols, _Opt> > >(mat) );
}

template <class T, int _Rows, int _Cols, int _Opt>
typename gsSolverOp< gsEigen::LDLT< gsEigen::Matrix<T, _Rows, _Cols, _Opt> > >::uPtr  makeCholeskySolver(const memory::shared_ptr< gsMatrix<T, _Rows, _Cols, _Opt> > & mat)
{
    return memory::make_unique( new gsSolverOp< gsEigen::LDLT< gsEigen::Matrix<T, _Rows, _Cols, _Opt> > >(mat) );
}


template <typename T, int _Opt, typename _Index>
typename gsSolverOp< typename gsSparseSolver<T>::LU >::uPtr  makeSparseLUSolver(const gsSparseMatrix<T,_Opt,_Index> & mat)
{
    return memory::make_unique( new gsSolverOp< typename gsSparseSolver<T>::LU >(mat) );
}

template <typename T, int _Opt, typename _Index>
typename gsSolverOp< typename gsSparseSolver<T>::LU >::uPtr  makeSparseLUSolver(const memory::shared_ptr< gsSparseMatrix<T,_Opt,_Index> > & mat)
{
    return memory::make_unique( new gsSolverOp< typename gsSparseSolver<T>::LU >(mat) );
}


template <typename T, int _Opt, typename _Index>
typename gsSolverOp< typename gsSparseSolver<T>::QR >::uPtr  makeSparseQRSolver(const gsSparseMatrix<T,_Opt,_Index> & mat)
{
    return memory::make_unique( new gsSolverOp< typename gsSparseSolver<T>::QR >(mat) );
}

template <typename T, int _Opt, typename _Index>
typename gsSolverOp< typename gsSparseSolver<T>::QR >::uPtr  makeSparseQRSolver(const memory::shared_ptr< gsSparseMatrix<T,_Opt,_Index> > & mat)
{
    return memory::make_unique( new gsSolverOp< typename gsSparseSolver<T>::QR >(mat) );
}

 
template <typename T, int _Opt, typename _Index>
typename gsSolverOp< typename gsSparseSolver<T>::SimplicialLDLT >::uPtr  makeSparseCholeskySolver(const gsSparseMatrix<T,_Opt,_Index> & mat)
{
    return memory::make_unique( new gsSolverOp<typename  gsSparseSolver<T>::SimplicialLDLT >(mat) );
}

template <typename T, int _Opt, typename _Index>
typename gsSolverOp< typename gsSparseSolver<T>::SimplicialLDLT >::uPtr  makeSparseCholeskySolver(const memory::shared_ptr< gsSparseMatrix<T,_Opt,_Index> > & mat)
{
    return memory::make_unique( new gsSolverOp<typename  gsSparseSolver<T>::SimplicialLDLT >(mat) );
}


} // namespace gismo