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>