Commit bac4cdfa authored by Tomáš Oberhuber's avatar Tomáš Oberhuber Committed by Tomáš Oberhuber
Browse files

Revision of multidiagonal matrix.

parent f40eb2d7

src/TNL/Matrices/Multidiagonal.h

0 → 100644
+217 −0
Original line number Diff line number Diff line 
/***************************************************************************

                          Multidiagonal.h  -  description

                             -------------------

    begin                : Oct 13, 2011

    copyright            : (C) 2011 by Tomas Oberhuber

    email                : tomas.oberhuber@fjfi.cvut.cz

 ***************************************************************************/



/* See Copyright Notice in tnl/Copyright */



#pragma once



#include <TNL/Matrices/Matrix.h>

#include <TNL/Containers/Vector.h>

#include <TNL/Matrices/MultidiagonalMatrixRowView.h>

#include <TNL/Containers/Segments/Ellpack.h>

#include <TNL/Matrices/details/MultidiagonalMatrixIndexer.h>

#include <TNL/Matrices/MultidiagonalMatrixView.h>



namespace TNL {

namespace Matrices {



template< typename Real = double,

          typename Device = Devices::Host,

          typename Index = int,

          bool RowMajorOrder = std::is_same< Device, Devices::Host >::value,

          typename RealAllocator = typename Allocators::Default< Device >::template Allocator< Real >,

          typename IndexAllocator = typename Allocators::Default< Device >::template Allocator< Index > >

class Multidiagonal : public Matrix< Real, Device, Index, RealAllocator >

{

   public:

      using RealType = Real;

      using DeviceType = Device;

      using IndexType = Index;

      using RealAllocatorType = RealAllocator;

      using IndexAllocatorType = IndexAllocator;

      using BaseType = Matrix< Real, Device, Index, RealAllocator >;

      using ValuesType = typename BaseType::ValuesVector;

      using ValuesViewType = typename ValuesType::ViewType;

      using IndexerType = details::MultidiagonalMatrixIndexer< IndexType, RowMajorOrder >;

      using RowView = MultidiagonalMatrixRowView< ValuesViewType, IndexerType >;

      using ViewType = MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >;

      using ConstViewType = MultidiagonalMatrixView< typename std::add_const< Real >::type, Device, Index, RowMajorOrder >;



      using DiagonalsShiftsType = Containers::Vector< IndexType, DeviceType, IndexType, IndexAllocatorType >;

      using DiagonalsShiftsView = typename DiagonalsShiftsType::ViewType;

      using HostDiagonalsShiftsType = Containers::Vector< IndexType, Devices::Host, IndexType >;

      using HostDiagonalsShiftsView = typename HostDiagonalsShiftsType::ViewType;





      // TODO: remove this - it is here only for compatibility with original matrix implementation

      typedef Containers::Vector< IndexType, DeviceType, IndexType > CompressedRowLengthsVector;

      typedef Containers::VectorView< IndexType, DeviceType, IndexType > CompressedRowLengthsVectorView;

      typedef typename CompressedRowLengthsVectorView::ConstViewType ConstCompressedRowLengthsVectorView;



      template< typename _Real = Real,

                typename _Device = Device,

                typename _Index = Index >

      using Self = Multidiagonal< _Real, _Device, _Index >;



      static constexpr bool getRowMajorOrder() { return RowMajorOrder; };



      Multidiagonal();



      Multidiagonal( const IndexType rows,

                     const IndexType columns );



      template< typename Vector >

      Multidiagonal( const IndexType rows,

                     const IndexType columns,

                     const Vector& diagonalsShifts );



      ViewType getView() const; // TODO: remove const



      //ConstViewType getConstView() const;



      static String getSerializationType();



      virtual String getSerializationTypeVirtual() const;



      template< typename Vector >

      void setDimensions( const IndexType rows,

                          const IndexType columns,

                          const Vector&  diagonalsShifts );



      //template< typename Vector >

      void setCompressedRowLengths( const ConstCompressedRowLengthsVectorView rowCapacities );



      const IndexType& getDiagonalsCount() const;



      const DiagonalsShiftsType& getDiagonalsShifts() const;



      template< typename Vector >

      void getCompressedRowLengths( Vector& rowLengths ) const;



      IndexType getNonemptyRowsCount() const;



      [[deprecated]]

      IndexType getRowLength( const IndexType row ) const;



      IndexType getMaxRowLength() const;



      template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_, typename RealAllocator_ >

      void setLike( const Multidiagonal< Real_, Device_, Index_, RowMajorOrder_, RealAllocator_ >& m );



      IndexType getNumberOfNonzeroMatrixElements() const;



      void reset();



      template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_, typename RealAllocator_ >

      bool operator == ( const Multidiagonal< Real_, Device_, Index_, RowMajorOrder_, RealAllocator_ >& matrix ) const;



      template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_, typename RealAllocator_ >

      bool operator != ( const Multidiagonal< Real_, Device_, Index_, RowMajorOrder_, RealAllocator_ >& matrix ) const;



      RowView getRow( const IndexType& rowIdx );



      const RowView getRow( const IndexType& rowIdx ) const;



      void setValue( const RealType& v );



      bool setElement( const IndexType row,

                       const IndexType column,

                       const RealType& value );



      bool addElement( const IndexType row,

                       const IndexType column,

                       const RealType& value,

                       const RealType& thisElementMultiplicator = 1.0 );



      RealType getElement( const IndexType row,

                           const IndexType column ) const;



      template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >

      void rowsReduction( IndexType first, IndexType last, Fetch& fetch, Reduce& reduce, Keep& keep, const FetchReal& zero ) const;



      template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >

      void allRowsReduction( Fetch& fetch, Reduce& reduce, Keep& keep, const FetchReal& zero ) const;



      template< typename Function >

      void forRows( IndexType first, IndexType last, Function& function ) const;



      template< typename Function >

      void forRows( IndexType first, IndexType last, Function& function );



      template< typename Function >

      void forAllRows( Function& function ) const;



      template< typename Function >

      void forAllRows( Function& function );



      template< typename Vector >

      __cuda_callable__

      typename Vector::RealType rowVectorProduct( const IndexType row,

                                                  const Vector& vector ) const;



      template< typename InVector,

                typename OutVector >

      void vectorProduct( const InVector& inVector,

                          OutVector& outVector ) const;



      template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_, typename RealAllocator_ >

      void addMatrix( const Multidiagonal< Real_, Device_, Index_, RowMajorOrder_, RealAllocator_ >& matrix,

                      const RealType& matrixMultiplicator = 1.0,

                      const RealType& thisMatrixMultiplicator = 1.0 );



      template< typename Real2, typename Index2 >

      void getTransposition( const Multidiagonal< Real2, Device, Index2 >& matrix,

                             const RealType& matrixMultiplicator = 1.0 );



      template< typename Vector1, typename Vector2 >

      __cuda_callable__

      void performSORIteration( const Vector1& b,

                                const IndexType row,

                                Vector2& x,

                                const RealType& omega = 1.0 ) const;



      // copy assignment

      Multidiagonal& operator=( const Multidiagonal& matrix );



      // cross-device copy assignment

      template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_, typename RealAllocator_ >

      Multidiagonal& operator=( const Multidiagonal< Real_, Device_, Index_, RowMajorOrder_, RealAllocator_ >& matrix );



      void save( File& file ) const;



      void load( File& file );



      void save( const String& fileName ) const;



      void load( const String& fileName );



      void print( std::ostream& str ) const;



      const IndexerType& getIndexer() const;



      IndexerType& getIndexer();



   protected:



      __cuda_callable__

      IndexType getElementIndex( const IndexType row,

                                 const IndexType localIdx ) const;



      DiagonalsShiftsType diagonalsShifts;



      HostDiagonalsShiftsType hostDiagonalsShifts;



      IndexerType indexer;



      ViewType view;

};



} // namespace Matrices

} // namespace TNL



#include <TNL/Matrices/Multidiagonal.hpp>

src/TNL/Matrices/Multidiagonal.hpp

0 → 100644
+909 −0

File added.

Preview size limit exceeded, changes collapsed.

src/TNL/Matrices/MultidiagonalMatrixRowView.h

0 → 100644
+59 −0
Original line number Diff line number Diff line 
/***************************************************************************

                          MultidiagonalMatrixRowView.h  -  description

                             -------------------

    begin                : Jan 11, 2020

    copyright            : (C) 2020 by oberhuber

    email                : tomas.oberhuber@fjfi.cvut.cz

 ***************************************************************************/



/* See Copyright Notice in tnl/Copyright */



#pragma once



namespace TNL {

namespace Matrices {   



template< typename ValuesView,

          typename Indexer >

class MultidiagonalMatrixRowView

{

   public:



      using RealType = typename ValuesView::RealType;

      using IndexType = typename ValuesView::IndexType;

      using ValuesViewType = ValuesView;

      using IndexerType = Indexer;



      __cuda_callable__

      MultidiagonalMatrixRowView( const IndexType rowIdx,

                                const ValuesViewType& values,

                                const IndexerType& indexer );



      __cuda_callable__

      IndexType getSize() const;



      __cuda_callable__

      const IndexType getColumnIndex( const IndexType localIdx ) const;



      __cuda_callable__

      const RealType& getValue( const IndexType localIdx ) const;



      __cuda_callable__

      RealType& getValue( const IndexType localIdx );



      __cuda_callable__

      void setElement( const IndexType localIdx,

                       const RealType& value );

   protected:



      IndexType rowIdx;



      ValuesViewType values;



      Indexer indexer;

};



} // namespace Matrices

} // namespace TNL



#include <TNL/Matrices/MultidiagonalMatrixRowView.hpp>

src/TNL/Matrices/MultidiagonalMatrixRowView.hpp

0 → 100644
+75 −0
Original line number Diff line number Diff line 
/***************************************************************************

                          MultidiagonalMatrixRowView.hpp  -  description

                             -------------------

    begin                : Jan 11, 2020

    copyright            : (C) 2020 by oberhuber

    email                : tomas.oberhuber@fjfi.cvut.cz

 ***************************************************************************/



/* See Copyright Notice in tnl/Copyright */



#pragma once



namespace TNL {

namespace Matrices {   



template< typename ValuesView, typename Indexer >

__cuda_callable__

MultidiagonalMatrixRowView< ValuesView, Indexer >::

MultidiagonalMatrixRowView( const IndexType rowIdx,

                          const ValuesViewType& values,

                          const IndexerType& indexer )

: rowIdx( rowIdx ), values( values ), indexer( indexer )

{

}



template< typename ValuesView, typename Indexer >

__cuda_callable__

auto

MultidiagonalMatrixRowView< ValuesView, Indexer >::

getSize() const -> IndexType

{

   return indexer.getRowSize();

}



template< typename ValuesView, typename Indexer >

__cuda_callable__

auto

MultidiagonalMatrixRowView< ValuesView, Indexer >::

getColumnIndex( const IndexType localIdx ) const -> const IndexType

{

   TNL_ASSERT_GE( localIdx, 0, "" );

   TNL_ASSERT_LT( localIdx, 3, "" );

   return rowIdx + localIdx - 1;

}



template< typename ValuesView, typename Indexer >

__cuda_callable__

auto

MultidiagonalMatrixRowView< ValuesView, Indexer >::

getValue( const IndexType localIdx ) const -> const RealType&

{

   return this->values[ this->indexer.getGlobalIndex( rowIdx, localIdx ) ];

}



template< typename ValuesView, typename Indexer >

__cuda_callable__

auto

MultidiagonalMatrixRowView< ValuesView, Indexer >::

getValue( const IndexType localIdx ) -> RealType&

{

   return this->values[ this->indexer.getGlobalIndex( rowIdx, localIdx ) ];

}



template< typename ValuesView, typename Indexer >

__cuda_callable__

void 

MultidiagonalMatrixRowView< ValuesView, Indexer >::

setElement( const IndexType localIdx,

            const RealType& value )

{

   this->values[ indexer.getGlobalIndex( rowIdx, localIdx ) ] = value;

}



} // namespace Matrices

} // namespace TNL

src/TNL/Matrices/MultidiagonalMatrixView.h

0 → 100644
+181 −0
Original line number Diff line number Diff line 
/***************************************************************************

                          MultidiagonalMatrixView.h  -  description

                             -------------------

    begin                : Jan 11, 2020

    copyright            : (C) 2020 by Tomas Oberhuber

    email                : tomas.oberhuber@fjfi.cvut.cz

 ***************************************************************************/



/* See Copyright Notice in tnl/Copyright */



#pragma once



#include <TNL/Matrices/MatrixView.h>

#include <TNL/Containers/Vector.h>

#include <TNL/Matrices/MultidiagonalMatrixRowView.h>

#include <TNL/Containers/Segments/Ellpack.h>

#include <TNL/Matrices/details/MultidiagonalMatrixIndexer.h>



namespace TNL {

namespace Matrices {



template< typename Real = double,

          typename Device = Devices::Host,

          typename Index = int,

          bool RowMajorOrder = std::is_same< Device, Devices::Host >::value >

class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >

{

   public:

      using RealType = Real;

      using DeviceType = Device;

      using IndexType = Index;

      using BaseType = MatrixView< Real, Device, Index >;

      using DiagonalsShiftsType = Containers::Vector< IndexType, DeviceType, IndexType >;

      using DiagonalsShiftsView = typename DiagonalsShiftsType::ViewType;

      using HostDiagonalsShiftsType = Containers::Vector< IndexType, Devices::Host, IndexType >;

      using HostDiagonalsShiftsView = typename DiagonalsShiftsType::ViewType;

      using IndexerType = details::MultidiagonalMatrixIndexer< IndexType, RowMajorOrder >;

      using ValuesViewType = typename BaseType::ValuesView;

      using ViewType = MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >;

      using ConstViewType = MultidiagonalMatrixView< typename std::add_const< Real >::type, Device, Index, RowMajorOrder >;

      using RowView = MultidiagonalMatrixRowView< ValuesViewType, IndexerType >;



      // TODO: remove this - it is here only for compatibility with original matrix implementation

      typedef Containers::Vector< IndexType, DeviceType, IndexType > CompressedRowLengthsVector;

      typedef Containers::VectorView< IndexType, DeviceType, IndexType > CompressedRowLengthsVectorView;

      typedef typename CompressedRowLengthsVectorView::ConstViewType ConstCompressedRowLengthsVectorView;



      template< typename _Real = Real,

                typename _Device = Device,

                typename _Index = Index,

                bool RowMajorOrder_ = std::is_same< Device, Devices::Host >::value >

      using Self = MultidiagonalMatrixView< _Real, _Device, _Index, RowMajorOrder_ >;



      MultidiagonalMatrixView();



      MultidiagonalMatrixView( const ValuesViewType& values,

                               const DiagonalsShiftsView& diagonalsShifts,

                               const IndexerType& indexer );



      ViewType getView();



      ConstViewType getConstView() const;



      static String getSerializationType();



      virtual String getSerializationTypeVirtual() const;



      __cuda_callable__

      const IndexType& getDiagonalsCount() const;



      template< typename Vector >

      void getCompressedRowLengths( Vector& rowLengths ) const;



      IndexType getNonemptyRowsCount() const;



      [[deprecated]]

      IndexType getRowLength( const IndexType row ) const;



      IndexType getMaxRowLength() const;



      IndexType getNumberOfNonzeroMatrixElements() const;



      template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_ >

      bool operator == ( const MultidiagonalMatrixView< Real_, Device_, Index_, RowMajorOrder_ >& matrix ) const;



      template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_ >

      bool operator != ( const MultidiagonalMatrixView< Real_, Device_, Index_, RowMajorOrder_ >& matrix ) const;



      RowView getRow( const IndexType& rowIdx );



      const RowView getRow( const IndexType& rowIdx ) const;



      void setValue( const RealType& v );



      bool setElement( const IndexType row,

                       const IndexType column,

                       const RealType& value );



      bool addElement( const IndexType row,

                       const IndexType column,

                       const RealType& value,

                       const RealType& thisElementMultiplicator = 1.0 );



      RealType getElement( const IndexType row,

                           const IndexType column ) const;



      MultidiagonalMatrixView& operator=( const MultidiagonalMatrixView& view );



      template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >

      void rowsReduction( IndexType first, IndexType last, Fetch& fetch, Reduce& reduce, Keep& keep, const FetchReal& zero ) const;



      template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >

      void allRowsReduction( Fetch& fetch, Reduce& reduce, Keep& keep, const FetchReal& zero ) const;



      template< typename Function >

      void forRows( IndexType first, IndexType last, Function& function ) const;



      template< typename Function >

      void forRows( IndexType first, IndexType last, Function& function );



      template< typename Function >

      void forAllRows( Function& function ) const;



      template< typename Function >

      void forAllRows( Function& function );



      template< typename Vector >

      __cuda_callable__

      typename Vector::RealType rowVectorProduct( const IndexType row,

                                                  const Vector& vector ) const;



      template< typename InVector,

                typename OutVector >

      void vectorProduct( const InVector& inVector,

                          OutVector& outVector ) const;



      template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_ >

      void addMatrix( const MultidiagonalMatrixView< Real_, Device_, Index_, RowMajorOrder_ >& matrix,

                      const RealType& matrixMultiplicator = 1.0,

                      const RealType& thisMatrixMultiplicator = 1.0 );



      template< typename Real2, typename Index2 >

      void getTransposition( const MultidiagonalMatrixView< Real2, Device, Index2 >& matrix,

                             const RealType& matrixMultiplicator = 1.0 );



      template< typename Vector1, typename Vector2 >

      __cuda_callable__

      void performSORIteration( const Vector1& b,

                                const IndexType row,

                                Vector2& x,

                                const RealType& omega = 1.0 ) const;



      void save( File& file ) const;



      void save( const String& fileName ) const;



      void print( std::ostream& str ) const;



      __cuda_callable__

      const IndexerType& getIndexer() const;



      __cuda_callable__

      IndexerType& getIndexer();



   protected:



      __cuda_callable__

      IndexType getElementIndex( const IndexType row,

                                 const IndexType localIdx ) const;



      DiagonalsShiftsView diagonalsShifts;



      HostDiagonalsShiftsView hostDiagonalsShifts;



      IndexerType indexer;

};



} // namespace Matrices

} // namespace TNL



#include <TNL/Matrices/MultidiagonalMatrixView.hpp>