gsAsMatrix.h

Go to the documentation of this file.

# pragma once

// Assumes that Eigen library has been already included

namespace gismo
{

template<class T, int _Rows, int _Cols>
class gsAsMatrix : public gsEigen::Map< gsEigen::Matrix<T,_Rows,_Cols> >
{
public:
    typedef gsEigen::Map< gsEigen::Matrix<T,_Rows,_Cols> > Base;

    // type of first minor matrix: rows and cols reduced by one
    typedef typename gsMatrix< T, _Rows, _Cols>::FirstMinorMatrixType FirstMinorMatrixType;
    
    // type of row minor matrix: rows reduced by one
    typedef typename gsMatrix< T, _Rows, _Cols>::RowMinorMatrixType RowMinorMatrixType;

    // type of col minor matrix: cols reduced by one
    typedef typename gsMatrix< T, _Rows, _Cols>::ColMinorMatrixType  ColMinorMatrixType;

    // Type pointing to the transpose of the matrix
    typedef gsEigen::Transpose<Base> Tr;

    // Type pointing to the (const) transpose of the matrix
    typedef const gsEigen::Transpose<const Base> constTr;

public:
    gsAsMatrix( std::vector<T> & v, index_t n, index_t m)
    : Base( v.data(), n, m)
    { 
        //GISMO_ASSERT( v.size() != 0, "Tried to map an empty vector." ); 
        GISMO_ASSERT( m*n <= (index_t)(v.size()), "Not enough coefficients in vector to map." );
    }

    gsAsMatrix( std::vector<T> & v)
    : Base( v.data(), 1, v.size() ) 
    {  
        //GISMO_ASSERT( v.size() != 0, "Tried to map an empty vector." ); 
    }

    gsAsMatrix( T * pt, unsigned n, unsigned m)
    : Base( pt, n, m) {  }

    gsMatrix<T> move()
    {
        gsMatrix<T> a;
        a.swap(*this);
        return a;
    }
    
#ifdef _MSC_VER
    template <class EigenExpr>
    gsAsMatrix& operator= (const EigenExpr & other) 
    {
        this->Base::operator=(other);
        return *this;
    }
#else
    using Base::operator=;
#endif

    void firstMinor(index_t i, index_t j, FirstMinorMatrixType & result ) const
    {
        const index_t mrows = this->rows()-1, 
            mcols = this->cols()-1;
        GISMO_ASSERT( i <= mrows, "Invalid row." );
        GISMO_ASSERT( j <= mcols, "Invalid column." );
        result.resize(mrows,mcols);
        result.block(0,0,i,j)             = this->block(0,0,i,j);
        result.block(i,0,mrows-i,j)       = this->block(i+1,0,mrows-i,j);
        result.block(0,j,i,mcols-j)       = this->block(0,j+1,i,mcols-j);
        result.block(i,j,mrows-i,mcols-j) = this->block(i+1,j+1,mrows-i,mcols-j);
    }

    void rowMinor(index_t i, RowMinorMatrixType & result ) const
    {
        const index_t mrows = this->rows()-1;
        GISMO_ASSERT( 0 <= i && i <= mrows, "Invalid row." );
        result.resize(mrows, this->cols());
        result.topRows(i)          = this->topRows(i);
        result.bottomRows(mrows-i) = this->bottomRows(mrows-i);
    }
    
    void colMinor(index_t j, ColMinorMatrixType & result ) const
    {
        const index_t mcols = this->cols()-1;
        GISMO_ASSERT( 0 <= j && j <= mcols, "Invalid column." );
        result.resize( this->rows(), mcols);
        result.leftCols(j)        = this->leftCols(j);
        result.rightCols(mcols-j) = this->rightCols(mcols-j);
    }

private:
    gsAsMatrix();
};

template<class T, int _Rows, int _Cols>
class gsAsConstMatrix : public gsEigen::Map< const gsEigen::Matrix<T,_Rows,_Cols> >
{
public:
    typedef gsEigen::Map<const gsEigen::Matrix<T,_Rows,_Cols> > Base;

    // Type pointing to the transpose of the matrix
    typedef gsEigen::Transpose<Base> Tr;

    // Type pointing to the (const) transpose of the matrix
    typedef const gsEigen::Transpose<const Base> constTr;

    // type of first minor matrix: rows and cols reduced by one
    typedef typename gsMatrix< T, _Rows, _Cols>::FirstMinorMatrixType FirstMinorMatrixType;
    
    // type of row minor matrix: rows reduced by one
    typedef typename gsMatrix< T, _Rows, _Cols>::RowMinorMatrixType RowMinorMatrixType;

    // type of col minor matrix: cols reduced by one
    typedef typename gsMatrix< T, _Rows, _Cols>::ColMinorMatrixType  ColMinorMatrixType;

public:

    gsAsConstMatrix( const std::vector<T> & v, index_t n, index_t m)
    : Base( v.data(), n, m)
    { 
        GISMO_ASSERT( m*n <= (index_t)(v.size()), "Not enough coefficients in vector to map." );
    }

    gsAsConstMatrix( const std::vector<T> & v)
    : Base( v.data(), 1, v.size() ) 
    {  
        //GISMO_ASSERT( v.size() != 0, "Tried to map an empty vector." ); 
    }

    gsAsConstMatrix( const T * pt, unsigned n, unsigned m)
    : Base( pt, n, m) {  }

    gsAsConstMatrix(const gsEigen::Map< gsEigen::Matrix<T,_Rows,_Cols> > & mat)
    : Base( mat.data(), mat.rows(), mat.cols()) 
    {  }

public:
    
    void firstMinor(index_t i, index_t j, FirstMinorMatrixType & result ) const
    {
        const index_t mrows = this->rows()-1, 
            mcols = this->cols()-1;
        GISMO_ASSERT( i <= mrows, "Invalid row." );
        GISMO_ASSERT( j <= mcols, "Invalid column." );
        result.resize(mrows,mcols);
        result.block(0,0,i,j)             = this->block(0,0,i,j);
        result.block(i,0,mrows-i,j)       = this->block(i+1,0,mrows-i,j);
        result.block(0,j,i,mcols-j)       = this->block(0,j+1,i,mcols-j);
        result.block(i,j,mrows-i,mcols-j) = this->block(i+1,j+1,mrows-i,mcols-j);
    }

    void rowMinor(index_t i, RowMinorMatrixType & result ) const
    {
        const index_t mrows = this->rows()-1;
        GISMO_ASSERT( 0 <= i && i <= mrows, "Invalid row." );
        result.resize(mrows, this->cols());
        result.topRows(i)          = this->topRows(i);
        result.bottomRows(mrows-i) = this->bottomRows(mrows-i);
    }
    
    void colMinor(index_t j, ColMinorMatrixType & result ) const
    {
        const index_t mcols = this->cols()-1;
        GISMO_ASSERT( 0 <= j && j <= mcols, "Invalid column." );
        result.resize( this->rows(), mcols);
        result.leftCols(j)        = this->leftCols(j);
        result.rightCols(mcols-j) = this->rightCols(mcols-j);
    }

private:
    gsAsConstMatrix() { }
};

template<class T, int _Rows>
class gsAsVector : public gsAsMatrix<T,_Rows,1>
//class gsAsVector : public gsAsMatrix<T,_Rows,(_Rows==1?1:0)>
{
public:
    //typedef gsEigen::Map< gsEigen::Matrix<T,_Rows,1> > Base;
    //typedef gsAsMatrix<T,_Rows,(_Rows==1?1:0)> Base;
    typedef gsAsMatrix<T,_Rows,1> Base;

    // Type for treating a vector as a permutation matrix
    typedef gsEigen::PermutationMatrix<_Rows> Permutation;

public:
    gsAsVector( std::vector<T> & v)
    : Base( v.data(), v.size(), 1 ) 
    {  
        //GISMO_ASSERT( v.size() != 0, "Tried to map an empty vector." ); 
    }

    gsAsVector( T * pt, unsigned n)
    : Base( pt, n, 1) {  }

#ifdef _MSC_VER
    template <class EigenExpr>
    gsAsVector& operator= (const EigenExpr & other) 
    {
        this->Base::operator=(other);
        return *this;
    }
#else
    using Base::operator=;
#endif

private:
    gsAsVector() { }
};

template<class T, int _Rows>
//class gsAsConstVector : public gsAsConstMatrix<T,_Rows,(_Rows==1?1:0)>
class gsAsConstVector : public gsAsConstMatrix<T,_Rows,1>
{
public:
    //typedef gsEigen::Map<const gsEigen::Matrix<T,_Rows,1> > Base;

    typedef gsAsConstMatrix<T,_Rows,1> Base;
    //typedef gsAsConstMatrix<T,_Rows,(_Rows==1?1:0)> Base;
public:

    gsAsConstVector( const std::vector<T> & v)
    : Base( v.data(), v.size(), 1) 
    {  
        //GISMO_ASSERT( v.size() != 0, "Tried to map an empty vector." ); 
    }

    gsAsConstVector( const T * pt, unsigned n)
    : Base( pt, n, 1) {  }

private:
    gsAsConstVector() { }
};


template<class T, class iterator>
gsMatrix<T> makeMatrix(iterator it, index_t n, index_t m)
{
    gsMatrix<T> result(n,m);
    for ( index_t i = 0; i!=n; ++i)
        for ( index_t j = 0; j!=m; ++j)
            result(i,j)= *(it++);
    return result;
}


} // namespace gismo