gsLinearAlgebra.h

Go to the documentation of this file.

# pragma once

#include <gsCore/gsMath.h>

// Eigen linear algebra library (http://eigen.tuxfamily.org)

#define Eigen gsEigen
// Make Eigen use GISMO_ASSERT which throws exceptions
//
// Must be defined before including Eigen headers
// http://eigen.tuxfamily.org/dox-3.2/TopicPreprocessorDirectives.html
#define eigen_assert( cond ) GISMO_ASSERT( cond, "" )

// Plugin provides extra members
#define EIGEN_MATRIXBASE_PLUGIN <gsMatrix/gsMatrixAddons.h>
#define EIGEN_PLAINOBJECTBASE_PLUGIN <gsMatrix/gsPlainObjectBaseAddons.h>
#include <gsMatrix/gsEigenDeclarations.h>

#include <gsEigen/Core>

#if defined(gsMpfr_ENABLED)
#include <unsupported/Eigen/MPRealSupport>
#endif

#if defined(gsGmp_ENABLED)
#include <unsupported/Eigen/MPQClassSupport>
#endif

#include <gsEigen/Dense>
#include <gsEigen/Sparse>
#include <gsEigen/StdVector>
#include <gsEigen/Geometry>

// Extra Eigen code
#include <gsMatrix/Adjugate.h>
#include <gsMatrix/BlockDiag.h>
#include <gsMatrix/BlockTranspose.h>
//#include <gsMatrix/RowSelection.h>

#ifdef GISMO_WITH_SUPERLU
#include <gsEigen/SuperLUSupport>
#endif

#ifdef GISMO_WITH_PARDISO
#include <gsEigen/PardisoSupport>
#endif

#ifdef GISMO_WITH_PASTIX
#include <gsEigen/PaStiXSupport>
#endif

// sparsesuite
//#include <gsEigen/UmfPackSupport>
//#include <gsEigen/SPQRSupport>
//#include <gsEigen/CholmodSupport>

// METIS
//#include <gsEigen/MetisSupport>

// PaStiX
//#include <gsEigen/PaStiXSupport>

#ifdef GISMO_WITH_PYBIND11
#include <pybind11/eigen.h>
#endif

#undef Eigen
#undef eigen_assert
#undef EIGEN_MATRIXBASE_PLUGIN
#undef EIGEN_PLAINOBJECTBASE_PLUGIN


namespace gismo
{

using gsEigen::internal::cast; // from Core/MathFunctions.h

template<typename Derived>
inline bool (isnumber)(const gsEigen::MatrixBase<Derived>& x)
{ return ((x.array() == x.array())).all(); }

template<typename Derived>
inline bool isfinite(const gsEigen::MatrixBase<Derived>& x)
{ return ( (x - x).array() == (x - x).array()).all(); }



//Constantss related to gsMatrix
//( see also external/Eigen/src/Core/util/Constants.h )
using gsEigen::Dynamic ;//=-1

using gsEigen::Lower;//=1
using gsEigen::Upper;//=2

// Values for matrix align options
using gsEigen::RowMajor;//=0
using gsEigen::ColMajor;//=1
using gsEigen::AutoAlign;//=0

template<class T, int _Rows, int _Cols> class gsAsMatrix ;
template<class T, int _Rows, int _Cols> class gsAsConstMatrix ;

template<class T, int _Rows> class gsAsVector ;
template<class T, int _Rows> class gsAsConstVector ;

// helper template for changing the dimension of a matrix
template <int Dim, int Change>
struct ChangeDim
{
    enum { D = Change+Dim<0 ? 0 : Dim + Change };
};
template <int Change>
struct ChangeDim<Dynamic, Change>
{
    enum { D = Dynamic };
};


template<typename T>
struct gsEigenAdaptor
{
public:
    // Note: IncompleteILU is not compatible with
    // gsEigen::ConjugateGradient because this preconditionner does not
    // preserve symmetry.

    typedef gsEigen::ConjugateGradient<gsEigen::SparseMatrix<T,0,index_t>,
            gsEigen::Lower|gsEigen::Upper, gsEigen::IdentityPreconditioner> CGIdentity;

    typedef gsEigen::ConjugateGradient<gsEigen::SparseMatrix<T,0,index_t>, 
            gsEigen::Lower|gsEigen::Upper, gsEigen::DiagonalPreconditioner<T> > CGDiagonal;

    typedef gsEigen::BiCGSTAB<gsEigen::SparseMatrix<T,0,index_t>,
                            gsEigen::IncompleteLUT<T, index_t> > BiCGSTABILUT;

    typedef gsEigen::BiCGSTAB<gsEigen::SparseMatrix<T,0,index_t>,
                            gsEigen::DiagonalPreconditioner<T> > BiCGSTABDiagonal;

    typedef gsEigen::BiCGSTAB<gsEigen::SparseMatrix<T,0,index_t>,
                            gsEigen::IdentityPreconditioner > BiCGSTABIdentity;

    typedef gsEigen::SimplicialLDLT<gsEigen::SparseMatrix<T,0,index_t> > SimplicialLDLT;

    typedef gsEigen::SimplicialLLT<gsEigen::SparseMatrix<T,0,index_t> > SimplicialLLT;

    typedef gsEigen::SparseLU<gsEigen::SparseMatrix<T,0,index_t>,
                            gsEigen::COLAMDOrdering<index_t> > SparseLU;

    typedef gsEigen::SparseQR<gsEigen::SparseMatrix<T,0,index_t>,
                            gsEigen::COLAMDOrdering<index_t> > SparseQR;
    
    #ifdef GISMO_WITH_SUPERLU
    typedef gsEigen::SuperLU<gsEigen::SparseMatrix<T,0,index_t> > SuperLU;
    #endif

    #ifdef GISMO_WITH_PARDISO
    typedef gsEigen::PardisoLDLT<gsEigen::SparseMatrix<T,0,int> > PardisoLDLT;
    typedef gsEigen::PardisoLLT <gsEigen::SparseMatrix<T,0,int> > PardisoLLT;
    typedef gsEigen::PardisoLU  <gsEigen::SparseMatrix<T,0,int> > PardisoLU;
    #endif

};

} // namespace gismo


#include <gsMatrix/gsMatrixBlockView.h>
#include <gsMatrix/gsMatrix.h>
#include <gsMatrix/gsVector.h>
#include <gsMatrix/gsAsMatrix.h>
#include <gsMatrix/gsSparseMatrix.h>
#include <gsMatrix/gsSparseVector.h>
#include <gsMatrix/gsSparseSolver.h>
#include <gsMatrix/gsPoint.h>