SparseBlockMatrixBase.hpp

/*
 * This file is part of bogus, a C++ sparse block matrix library.
 *
 * Copyright 2013 Gilles Daviet <gdaviet@gmail.com>
 *
 * This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, You can obtain one at http://mozilla.org/MPL/2.0/.
*/


#ifndef BOGUS_SPARSEBLOCKMATRIX_BASE_HPP
#define BOGUS_SPARSEBLOCKMATRIX_BASE_HPP

#include "BlockMatrixBase.hpp"
#include "Expressions.hpp"
#include "DynamicExpressions.hpp"

#include "SparseBlockIndex.hpp"
#include "CompressedSparseBlockIndex.hpp"

#include "../Utils/CppTools.hpp"
#include "../Utils/Lock.hpp"

namespace bogus
{

template < typename BlockT, int Flags > class SparseBlockMatrix  ;
template < bool Symmetric > struct SparseBlockMatrixFinalizer ;
template < typename Derived, bool Major > struct SparseBlockIndexGetter ;


template < typename Derived >
class SparseBlockMatrixBase : public BlockMatrixBase< Derived >
{

public:
    // Convenient typedefs  and using directives

    typedef BlockMatrixBase< Derived > Base ;
    typedef BlockMatrixTraits< Derived > Traits ;

    typedef typename Traits::BlockPtr               BlockPtr ;
    typedef typename Base::Index                    Index ;


    typedef typename Traits::MajorIndexType         MajorIndexType ;

    typedef typename Base::BlockType                BlockType ;
    typedef typename Base::Scalar                   Scalar ;

    typedef typename MajorIndexType::InnerIterator  InnerIterator ;

    typedef SparseBlockIndex< false, Index, BlockPtr > UncompressedIndexType ;
    typedef SparseBlockIndex<  true, Index, BlockPtr > CompressedIndexType ;

    // Minor index is always uncompressed, as the blocks cannot be contiguous
    // For a symmetric matrix, it does not store diagonal block in the minor and transpose index
    typedef UncompressedIndexType MinorIndexType ;
    // Transpose index is compressed for perf, as we always create it in a compressed-compatible way
    typedef CompressedIndexType TransposeIndexType ;

    typedef typename TypeSwapIf< Traits::is_col_major, MajorIndexType, MinorIndexType >::First  RowIndexType ;
    typedef typename TypeSwapIf< Traits::is_col_major, MajorIndexType, MinorIndexType >::Second ColIndexType ;

    // Canonical type for a mutable matrix with different block type
    template < typename OtherBlockType, bool PreserveSymmetry = true, bool SwitchDirection = false >
    struct MutableImpl
    {
        typedef SparseBlockMatrix< OtherBlockType,
            Traits::flags & ~flags::UNCOMPRESSED > Type ;
    } ;
    template < typename OtherBlockType >
    struct MutableImpl< OtherBlockType, false, false >
    {
        typedef SparseBlockMatrix< OtherBlockType,
            Traits::flags & ~flags::UNCOMPRESSED & ~flags::SYMMETRIC > Type ;
    } ;
    template < typename OtherBlockType >
    struct MutableImpl< OtherBlockType, true, true >
    {
        typedef SparseBlockMatrix< OtherBlockType,
            (Traits::flags & ~flags::UNCOMPRESSED &
            ~flags::COL_MAJOR) | ((~Traits::flags)&flags::COL_MAJOR) > Type ;
    } ;
    template < typename OtherBlockType >
    struct MutableImpl< OtherBlockType, false, true >
    {
        typedef SparseBlockMatrix< OtherBlockType,
            (Traits::flags & ~flags::UNCOMPRESSED & ~flags::SYMMETRIC&
            ~flags::COL_MAJOR) | ((~Traits::flags)&flags::COL_MAJOR) > Type ;
    } ;
    typedef typename MutableImpl< BlockType, true >::Type CopyResultType ;

    typedef typename Base::ConstTransposeReturnType ConstTransposeReturnType ;

    typedef typename BlockTraits< BlockType >::TransposeStorageType TransposeBlockType ;
    typedef typename ResizableSequenceContainer< TransposeBlockType >::Type TransposeArrayType ;

    using Base::rows ;
    using Base::cols ;
    using Base::blocks ;
    using Base::derived ;
    using Base::block ;
    using Base::diagonal ;
    using Base::InvalidBlockPtr ;

    // Public API



    void setRows( const Index nBlocks, const unsigned* rowsPerBlock ) ;
    void setRows( const std::vector< unsigned > &rowsPerBlock )
    {
        setRows( rowsPerBlock.size(), &rowsPerBlock[0] ) ;
    }
    void setRows( const Index nBlocks, const Index rowsPerBlock )
    {
        setRows( std::vector< unsigned >( nBlocks, rowsPerBlock ) ) ;
    }
    void setRows( const Index nBlocks )
    {
        assert( Base::has_fixed_rows_blocks ) ;
        setRows( nBlocks, BlockType::RowsAtCompileTime ) ;
    }


    void setCols( const Index nBlocks, const unsigned* colsPerBlock ) ;
    void setCols( const std::vector< unsigned > &colsPerBlock )
    {
        setCols( colsPerBlock.size(), &colsPerBlock[0] ) ;
    }
    void setCols( const Index nBlocks, const Index colsPerBlock )
    {
        setCols( std::vector< unsigned >( nBlocks, colsPerBlock ) ) ;
    }
    void setCols( const Index nBlocks )
    {
        assert( Base::has_fixed_cols_blocks ) ;
        setCols( nBlocks, BlockType::ColsAtCompileTime ) ;
    }

    Index rowsOfBlocks() const { return rowMajorIndex().outerSize() ; }
    Index colsOfBlocks() const { return colMajorIndex().outerSize() ; }

    Index blockRows( Index row ) const { return rowOffsets()[ row + 1 ] - rowOffsets()[ row ] ; }
    Index blockCols( Index col ) const { return colOffsets()[ col + 1 ] - colOffsets()[ col ] ; }



    void reserve( std::size_t nBlocks )
    {
        m_blocks.reserve( nBlocks ) ;
        m_majorIndex.reserve( nBlocks ) ;
    }


    BlockType& insertBack( Index row, Index col )
    {
        if( Traits::is_col_major )
            return insertByOuterInner< true >( col, row ) ;
        else
            return insertByOuterInner< true >( row, col ) ;
    }

    BlockType& insertBackAndResize( Index row, Index col ) ;


    BlockType& insert( Index row, Index col )
    {
        BOGUS_STATIC_ASSERT( !Traits::is_compressed, UNORDERED_INSERTION_WITH_COMPRESSED_INDEX ) ;
        if( Traits::is_col_major )
            return insertByOuterInner< false >( col, row ) ;
        else
            return insertByOuterInner< false >( row, col ) ;
    }

    BlockType& insertAndResize( Index row, Index col ) ;


    template< bool Ordered >
    BlockType& insertByOuterInner( Index outer, Index inner ) ;

    void finalize() ;

    void clear() ;
    Derived& setZero() { clear() ; return derived() ; }
    Derived& setIdentity() ;
    Derived& setBlocksToZero() ;

    std::size_t nBlocks() const { return blocks().size() ; }
    std::size_t size() const { return nBlocks() ; }

    bool empty() const { return 0 == nBlocks() ; }

    const BlockType& block( BlockPtr ptr ) const
    {
        return m_blocks[ ptr ] ;
    }

    BlockType& block( BlockPtr ptr )
    {
        return m_blocks[ ptr ] ;
    }

    BlockPtr blockPtr( Index row, Index col ) const ;
    BlockPtr diagonalBlockPtr( Index row  ) const ;

    template <bool ColWise, typename Func>
    void eachBlockOf( const Index outer, Func func ) const ;




    bool computeMinorIndex() ;


    void cacheTranspose() ;
    bool transposeCached() const { return m_transposeIndex.valid ; }


    InnerIterator innerIterator( Index outer ) const
    {
        return InnerIterator( m_majorIndex, outer ) ;
    }

    const MajorIndexType& majorIndex() const
    {
        return m_majorIndex ;
    }
    const MinorIndexType& minorIndex() const
    {
        return m_minorIndex ;
    }
    const TransposeIndexType& transposeIndex() const
    {
        return m_transposeIndex ;
    }
    const ColIndexType& colMajorIndex() const ;
    const RowIndexType& rowMajorIndex() const ;



    template < bool Transpose, typename OtherDerived >
    Derived& assign ( const SparseBlockMatrixBase< OtherDerived > &source, const Scalar scale = 1 ) ;

    template < typename OtherDerived >
    Derived& operator= ( const BlockObjectBase< OtherDerived > &source )
    {
        Evaluator< OtherDerived > rhs ( source.derived() ) ;
        return assign< OtherDerived::is_transposed >( *rhs ) ;
    }

    Derived& operator= ( const SparseBlockMatrixBase &source ) ;

    template < typename LhsT, typename RhsT >
    Derived& operator= ( const Product< LhsT, RhsT > &prod ) ;

    template < typename LhsT, typename RhsT >
    Derived& operator= ( const Addition< LhsT, RhsT > &add ) ;

    template < typename Expression >
    Derived& operator= ( const NarySum< Expression > &sum ) ;

    template < typename OtherDerived >
    Derived& operator= ( const Scaling< OtherDerived > &scaling )
    {
        Evaluator< typename Scaling< OtherDerived >::Operand::ObjectType > rhs ( scaling.operand.object ) ;
        return assign< Scaling< OtherDerived >::transposeOperand >( *rhs, scaling.operand.scaling ) ;
    }

    template< typename OtherDerived >
    void cloneDimensions( const BlockMatrixBase< OtherDerived > &source ) ;

    template< typename OtherDerived >
    void cloneStructure( const SparseBlockMatrixBase< OtherDerived > &source ) ;

    void cloneIndex( const MajorIndexType &index ) ;




    ConstTransposeReturnType transpose() const { return Transpose< Derived >( derived() ) ; }

    template < bool Transpose, typename RhsT, typename ResT >
    void multiply( const RhsT& rhs, ResT& res, Scalar alpha = 1, Scalar beta = 0 ) const ;

    template < typename RhsT, typename ResT >
    void splitRowMultiply( const Index row, const RhsT& rhs, ResT& res ) const ;

    template < bool DoTranspose, typename RhsT, typename ResT, typename PreOp >
    void rowMultiplyPrecompose( const Index row, const RhsT& rhs, ResT& res, const PreOp &op ) const ;

    template < bool DoTranspose, typename RhsT, typename ResT, typename PostOp >
    void colMultiplyPostcompose( const Index row, const RhsT& rhs, ResT& res, const PostOp &op ) const ;

    template < bool ColWise, typename LhsT, typename RhsT >
    void setFromProduct( const Product< LhsT, RhsT > &prod ) ;

    Derived& scale( Scalar alpha ) ;

    template < bool Transpose, typename OtherDerived >
    Derived& add( const SparseBlockMatrixBase< OtherDerived > &rhs, Scalar alpha = 1) ;

    Derived& operator *= ( Scalar alpha ) { return scale( alpha ) ; }
    Derived& operator /= ( Scalar alpha ) { return scale( 1./alpha ) ; }

    Scaling< Derived > operator-() const
    { return Scaling< Derived >( derived(), -1 ) ; }

    template < typename OtherDerived >
    Derived& operator+= ( const BlockObjectBase< OtherDerived > &source )
    {
        Evaluator< OtherDerived > rhs ( source.derived() ) ;
        return add< OtherDerived::is_transposed >( *rhs ) ;
    }

    template < typename OtherDerived >
    Derived& operator-= ( const BlockObjectBase< OtherDerived > &source )
    {
        Evaluator< OtherDerived > rhs ( source.derived() ) ;
        return add< OtherDerived::is_transposed >( *rhs, -1 ) ;
    }




#ifdef BOGUS_WITH_BOOST_SERIALIZATION
    template < typename Archive >
    void serialize( Archive & ar, const unsigned int file_version ) ;
#endif




    Derived& prune( const Scalar precision = 0 ) ;


    Derived& applyPermutation( const std::size_t* indices ) ;



    CopyResultType Zero() const
    {
        CopyResultType m ;
        m.cloneDimensions( *this ) ;
        return m.setZero() ;
    }

    CopyResultType Identity() const
    {
        CopyResultType m ;
        m.cloneDimensions( *this ) ;
        return m.setIdentity() ;
    }




    MajorIndexType& majorIndex() { return m_majorIndex ; }

    void prealloc( std::size_t nBlocks ) ;

    const TransposeArrayType& transposeBlocks() const
    { return m_transposeBlocks ; }

    const TransposeBlockType* transposeData() const
    { return data_pointer(m_transposeBlocks) ; }

    const Index* rowOffsets() const { return colMajorIndex().innerOffsetsData() ; }

    const Index* colOffsets() const { return rowMajorIndex().innerOffsetsData() ; }

    const Lock& lock() const { return m_lock ; }


protected:

    enum {
        is_bsr_compatible =
        Base::has_fixed_size_blocks && Base::has_square_or_dynamic_blocks &&
        Traits::is_compressed && !Traits::is_col_major &&
        BlockTraits< BlockType >::uses_plain_array_storage
    } ;

    using Base::m_cols ;
    using Base::m_rows ;
    using Base::m_blocks ;

    typedef SparseBlockMatrixFinalizer< Traits::is_symmetric > Finalizer ;
    friend struct SparseBlockIndexGetter< Derived, true > ;
    friend struct SparseBlockIndexGetter< Derived, false > ;

    SparseBlockMatrixBase() ;

    template< bool EnforceThreadSafety >
    BlockType& allocateBlock( BlockPtr &ptr ) ;

    void computeMinorIndex( MinorIndexType &cmIndex) const ;

    const MinorIndexType& getOrComputeMinorIndex( MinorIndexType &tempIndex) const ;

    ColIndexType& colMajorIndex() ;
    RowIndexType& rowMajorIndex() ;

    template< typename IndexT >
    void setInnerOffets( IndexT& index, const Index nBlocks, const unsigned* blockSizes ) const ;

    TransposeBlockType* transposeData()
    { return data_pointer(m_transposeBlocks) ; }


    TransposeArrayType m_transposeBlocks ;

    MajorIndexType m_majorIndex ;

    MinorIndexType m_minorIndex ;

    TransposeIndexType   m_transposeIndex ;

    Lock m_lock ;
} ;

}

#endif // SPARSEBLOCKMATRIX_HH