diff --git a/src/TNL/Matrices/Multidiagonal.h b/src/TNL/Matrices/Multidiagonal.h
new file mode 100644
index 0000000000000000000000000000000000000000..5d23cd960e5c69bdf718ff4e0a609757742c2526
--- /dev/null
+++ b/src/TNL/Matrices/Multidiagonal.h
@@ -0,0 +1,217 @@
+/***************************************************************************
+ 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>
diff --git a/src/TNL/Matrices/Multidiagonal.hpp b/src/TNL/Matrices/Multidiagonal.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..95f6667c1c1e99265d4b9714ff43f28807bc48f9
--- /dev/null
+++ b/src/TNL/Matrices/Multidiagonal.hpp
@@ -0,0 +1,909 @@
+/***************************************************************************
+ Multidiagonal.hpp - 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 <sstream>
+#include <TNL/Assert.h>
+#include <TNL/Matrices/Multidiagonal.h>
+#include <TNL/Exceptions/NotImplementedError.h>
+
+namespace TNL {
+namespace Matrices {
+
+template< typename Device >
+class MultidiagonalDeviceDependentCode;
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+Multidiagonal()
+{
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+ template< typename Vector >
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+Multidiagonal( const IndexType rows,
+ const IndexType columns,
+ const Vector& diagonalsShifts )
+{
+ this->setDimensions( rows, columns, diagonalsShifts );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+auto
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+getView() const -> ViewType
+{
+ // TODO: fix when getConstView works
+ return ViewType( const_cast< Multidiagonal* >( this )->values.getView(),
+ const_cast< Multidiagonal* >( this )->diagonalsShifts.getView(),
+ indexer );
+}
+
+/*template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+auto
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+getConstView() const -> ConstViewType
+{
+ return ConstViewType( this->values.getConstView(), indexer );
+}*/
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+String
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+getSerializationType()
+{
+ return String( "Matrices::Multidiagonal< " ) +
+ TNL::getSerializationType< RealType >() + ", [any_device], " +
+ TNL::getSerializationType< IndexType >() + ", " +
+ ( RowMajorOrder ? "true" : "false" ) + ", [any_allocator], [any_allocator] >";
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+String
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+getSerializationTypeVirtual() const
+{
+ return this->getSerializationType();
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+ template< typename Vector >
+void
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+setDimensions( const IndexType rows,
+ const IndexType columns,
+ const Vector& diagonalsShifts )
+{
+ Matrix< Real, Device, Index >::setDimensions( rows, columns );
+ this->diagonalsShifts = diagonalsShifts;
+ this->hostDiagonalsShifts = diagonalsShifts;
+ const IndexType minShift = min( diagonalsShifts );
+ IndexType nonemptyRows = min( rows, columns );
+ if( rows > columns && minShift < 0 )
+ nonemptyRows = min( rows, nonemptyRows - minShift );
+ this->indexer.set( rows, columns, diagonalsShifts.getSize(), nonemptyRows );
+ this->values.setSize( this->indexer.getStorageSize() );
+ this->values = 0.0;
+ this->view = this->getView();
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+ // template< typename Vector >
+void
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+setCompressedRowLengths( const ConstCompressedRowLengthsVectorView rowLengths )
+{
+ if( max( rowLengths ) > 3 )
+ throw std::logic_error( "Too many non-zero elements per row in a tri-diagonal matrix." );
+ if( rowLengths.getElement( 0 ) > 2 )
+ throw std::logic_error( "Too many non-zero elements per row in a tri-diagonal matrix." );
+ const IndexType diagonalLength = min( this->getRows(), this->getColumns() );
+ if( this->getRows() > this->getColumns() )
+ if( rowLengths.getElement( this->getRows()-1 ) > 1 )
+ throw std::logic_error( "Too many non-zero elements per row in a tri-diagonal matrix." );
+ if( this->getRows() == this->getColumns() )
+ if( rowLengths.getElement( this->getRows()-1 ) > 2 )
+ throw std::logic_error( "Too many non-zero elements per row in a tri-diagonal matrix." );
+ if( this->getRows() < this->getColumns() )
+ if( rowLengths.getElement( this->getRows()-1 ) > 3 )
+ throw std::logic_error( "Too many non-zero elements per row in a tri-diagonal matrix." );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+const Index&
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+getDiagonalsCount() const
+{
+ return this->view.getDiagonalsCount();
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+auto
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+getDiagonalsShifts() const -> const DiagonalsShiftsType&
+{
+ return this->diagonalsShifts;
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+ template< typename Vector >
+void
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+getCompressedRowLengths( Vector& rowLengths ) const
+{
+ return this->view.getCompressedRowLengths( rowLengths );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+Index
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+getNonemptyRowsCount() const
+{
+ return this->indexer.getNonemptyRowsCount();
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+Index
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+getRowLength( const IndexType row ) const
+{
+ return this->view.getRowLength( row );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+Index
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+getMaxRowLength() const
+{
+ return this->view.getMaxRowLength();
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+ template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_, typename RealAllocator_ >
+void
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+setLike( const Multidiagonal< Real_, Device_, Index_, RowMajorOrder_, RealAllocator_ >& m )
+{
+ this->setDimensions( m.getRows(), m.getColumns(), m.getDiagonalsShifts() );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+Index
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+getNumberOfNonzeroMatrixElements() const
+{
+ return this->view.getNumberOfNonzeroMatrixElements();
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+void
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+reset()
+{
+ Matrix< Real, Device, Index >::reset();
+ this->values.reset();
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+ template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_, typename RealAllocator_ >
+bool
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+operator == ( const Multidiagonal< Real_, Device_, Index_, RowMajorOrder_, RealAllocator_ >& matrix ) const
+{
+ if( RowMajorOrder == RowMajorOrder_ )
+ return this->values == matrix.values;
+ else
+ {
+ TNL_ASSERT( false, "TODO" );
+ }
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+ template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_, typename RealAllocator_ >
+bool
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+operator != ( const Multidiagonal< Real_, Device_, Index_, RowMajorOrder_, RealAllocator_ >& matrix ) const
+{
+ return ! this->operator==( matrix );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+void
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+setValue( const RealType& v )
+{
+ this->view.setValue( v );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+__cuda_callable__
+auto
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+getRow( const IndexType& rowIdx ) const -> const RowView
+{
+ return this->view.getRow( rowIdx );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+__cuda_callable__
+auto
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+getRow( const IndexType& rowIdx ) -> RowView
+{
+ return this->view.getRow( rowIdx );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+bool
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+setElement( const IndexType row, const IndexType column, const RealType& value )
+{
+ return this->view.setElement( row, column, value );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+bool
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+addElement( const IndexType row,
+ const IndexType column,
+ const RealType& value,
+ const RealType& thisElementMultiplicator )
+{
+ return this->view.addElement( row, column, value, thisElementMultiplicator );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+Real
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+getElement( const IndexType row, const IndexType column ) const
+{
+ return this->view.getElement( row, column );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+ template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >
+void
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+rowsReduction( IndexType first, IndexType last, Fetch& fetch, Reduce& reduce, Keep& keep, const FetchReal& zero ) const
+{
+ this->view.rowsReduction( first, last, fetch, reduce, keep, zero );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+ template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >
+void
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+allRowsReduction( Fetch& fetch, Reduce& reduce, Keep& keep, const FetchReal& zero ) const
+{
+ this->view.rowsReduction( 0, this->getRows(), fetch, reduce, keep, zero );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+ template< typename Function >
+void
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+forRows( IndexType first, IndexType last, Function& function ) const
+{
+ this->view.forRows( first, last, function );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+ template< typename Function >
+void
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+forRows( IndexType first, IndexType last, Function& function )
+{
+ this->view.forRows( first, last, function );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+ template< typename Function >
+void
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+forAllRows( Function& function ) const
+{
+ this->view.forRows( 0, this->getRows(), function );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+ template< typename Function >
+void
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+forAllRows( Function& function )
+{
+ this->view.forRows( 0, this->getRows(), function );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+template< typename Vector >
+__cuda_callable__
+typename Vector::RealType
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+rowVectorProduct( const IndexType row, const Vector& vector ) const
+{
+ return this->view.rowVectorProduct();
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+ template< typename InVector,
+ typename OutVector >
+void
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+vectorProduct( const InVector& inVector, OutVector& outVector ) const
+{
+ this->view.vectorProduct( inVector, outVector );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+ template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_, typename RealAllocator_ >
+void
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+addMatrix( const Multidiagonal< Real_, Device_, Index_, RowMajorOrder_, RealAllocator_ >& matrix,
+ const RealType& matrixMultiplicator,
+ const RealType& thisMatrixMultiplicator )
+{
+ this->view.addMatrix( matrix.getView(), matrixMultiplicator, thisMatrixMultiplicator );
+}
+
+#ifdef HAVE_CUDA
+template< typename Real,
+ typename Real2,
+ typename Index,
+ typename Index2 >
+__global__ void MultidiagonalTranspositionCudaKernel( const Multidiagonal< Real2, Devices::Cuda, Index2 >* inMatrix,
+ Multidiagonal< Real, Devices::Cuda, Index >* outMatrix,
+ const Real matrixMultiplicator,
+ const Index gridIdx )
+{
+ const Index rowIdx = ( gridIdx * Cuda::getMaxGridSize() + blockIdx.x ) * blockDim.x + threadIdx.x;
+ if( rowIdx < inMatrix->getRows() )
+ {
+ if( rowIdx > 0 )
+ outMatrix->setElementFast( rowIdx-1,
+ rowIdx,
+ matrixMultiplicator * inMatrix->getElementFast( rowIdx, rowIdx-1 ) );
+ outMatrix->setElementFast( rowIdx,
+ rowIdx,
+ matrixMultiplicator * inMatrix->getElementFast( rowIdx, rowIdx ) );
+ if( rowIdx < inMatrix->getRows()-1 )
+ outMatrix->setElementFast( rowIdx+1,
+ rowIdx,
+ matrixMultiplicator * inMatrix->getElementFast( rowIdx, rowIdx+1 ) );
+ }
+}
+#endif
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+ template< typename Real2, typename Index2 >
+void Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::getTransposition( const Multidiagonal< Real2, Device, Index2 >& matrix,
+ const RealType& matrixMultiplicator )
+{
+ TNL_ASSERT( this->getRows() == matrix.getRows(),
+ std::cerr << "This matrix rows: " << this->getRows() << std::endl
+ << "That matrix rows: " << matrix.getRows() << std::endl );
+ if( std::is_same< Device, Devices::Host >::value )
+ {
+ const IndexType& rows = matrix.getRows();
+ for( IndexType i = 1; i < rows; i++ )
+ {
+ RealType aux = matrix. getElement( i, i - 1 );
+ this->setElement( i, i - 1, matrix.getElement( i - 1, i ) );
+ this->setElement( i, i, matrix.getElement( i, i ) );
+ this->setElement( i - 1, i, aux );
+ }
+ }
+ if( std::is_same< Device, Devices::Cuda >::value )
+ {
+#ifdef HAVE_CUDA
+ Multidiagonal* kernel_this = Cuda::passToDevice( *this );
+ typedef Multidiagonal< Real2, Device, Index2 > InMatrixType;
+ InMatrixType* kernel_inMatrix = Cuda::passToDevice( matrix );
+ dim3 cudaBlockSize( 256 ), cudaGridSize( Cuda::getMaxGridSize() );
+ const IndexType cudaBlocks = roundUpDivision( matrix.getRows(), cudaBlockSize.x );
+ const IndexType cudaGrids = roundUpDivision( cudaBlocks, Cuda::getMaxGridSize() );
+ for( IndexType gridIdx = 0; gridIdx < cudaGrids; gridIdx++ )
+ {
+ if( gridIdx == cudaGrids - 1 )
+ cudaGridSize.x = cudaBlocks % Cuda::getMaxGridSize();
+ MultidiagonalTranspositionCudaKernel<<< cudaGridSize, cudaBlockSize >>>
+ ( kernel_inMatrix,
+ kernel_this,
+ matrixMultiplicator,
+ gridIdx );
+ }
+ Cuda::freeFromDevice( kernel_this );
+ Cuda::freeFromDevice( kernel_inMatrix );
+ TNL_CHECK_CUDA_DEVICE;
+#endif
+ }
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+ template< typename Vector1, typename Vector2 >
+__cuda_callable__
+void Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::performSORIteration( const Vector1& b,
+ const IndexType row,
+ Vector2& x,
+ const RealType& omega ) const
+{
+ RealType sum( 0.0 );
+ if( row > 0 )
+ sum += this->getElementFast( row, row - 1 ) * x[ row - 1 ];
+ if( row < this->getColumns() - 1 )
+ sum += this->getElementFast( row, row + 1 ) * x[ row + 1 ];
+ x[ row ] = ( 1.0 - omega ) * x[ row ] + omega / this->getElementFast( row, row ) * ( b[ row ] - sum );
+}
+
+
+// copy assignment
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >&
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::operator=( const Multidiagonal& matrix )
+{
+ this->setLike( matrix );
+ this->values = matrix.values;
+ return *this;
+}
+
+// cross-device copy assignment
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+ template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_, typename RealAllocator_ >
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >&
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+operator=( const Multidiagonal< Real_, Device_, Index_, RowMajorOrder_, RealAllocator_ >& matrix )
+{
+ static_assert( std::is_same< Device, Devices::Host >::value || std::is_same< Device, Devices::Cuda >::value,
+ "unknown device" );
+ static_assert( std::is_same< Device_, Devices::Host >::value || std::is_same< Device_, Devices::Cuda >::value,
+ "unknown device" );
+
+ this->setLike( matrix );
+ if( RowMajorOrder == RowMajorOrder_ )
+ this->values = matrix.getValues();
+ else
+ {
+ if( std::is_same< Device, Device_ >::value )
+ {
+ const auto matrix_view = matrix.getView();
+ auto f = [=] __cuda_callable__ ( const IndexType& rowIdx, const IndexType& localIdx, const IndexType& column, Real& value ) mutable {
+ value = matrix_view.getValues()[ matrix_view.getIndexer().getGlobalIndex( rowIdx, localIdx ) ];
+ };
+ this->forAllRows( f );
+ }
+ else
+ {
+ Multidiagonal< Real, Device, Index, RowMajorOrder_ > auxMatrix;
+ auxMatrix = matrix;
+ const auto matrix_view = auxMatrix.getView();
+ auto f = [=] __cuda_callable__ ( const IndexType& rowIdx, const IndexType& localIdx, const IndexType& column, Real& value ) mutable {
+ value = matrix_view.getValues()[ matrix_view.getIndexer().getGlobalIndex( rowIdx, localIdx ) ];
+ };
+ this->forAllRows( f );
+ }
+ }
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+void Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::save( File& file ) const
+{
+ Matrix< Real, Device, Index >::save( file );
+ file << diagonalsShifts;
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+void Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::load( File& file )
+{
+ Matrix< Real, Device, Index >::load( file );
+ file >> this->diagonalsShifts;
+ this->hostDiagonalsShifts = this->diagonalsShifts;
+ const IndexType minShift = min( diagonalsShifts );
+ IndexType nonemptyRows = min( this->getRows(), this->getColumns() );
+ if( this->getRows() > this->getColumns() && minShift < 0 )
+ nonemptyRows = min( this->getRows(), nonemptyRows - minShift );
+ this->indexer.set( this->getRows(), this->getColumns(), diagonalsShifts.getSize(), nonemptyRows );
+ this->view = this->getView();
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+void Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::save( const String& fileName ) const
+{
+ Object::save( fileName );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+void Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::load( const String& fileName )
+{
+ Object::load( fileName );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+void
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+print( std::ostream& str ) const
+{
+ for( IndexType row = 0; row < this->getRows(); row++ )
+ {
+ str <<"Row: " << row << " -> ";
+ for( IndexType column = row - 1; column < row + 2; column++ )
+ if( column >= 0 && column < this->columns )
+ str << " Col:" << column << "->" << this->getElement( row, column ) << "\t";
+ str << std::endl;
+ }
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+auto
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+getIndexer() const -> const IndexerType&
+{
+ return this->indexer;
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+auto
+Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+getIndexer() -> IndexerType&
+{
+ return this->indexer;
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder,
+ typename RealAllocator,
+ typename IndexAllocator >
+__cuda_callable__
+Index Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
+getElementIndex( const IndexType row, const IndexType column ) const
+{
+ IndexType localIdx = column - row;
+ if( row > 0 )
+ localIdx++;
+
+ TNL_ASSERT_GE( localIdx, 0, "" );
+ TNL_ASSERT_LT( localIdx, 3, "" );
+
+ return this->indexer.getGlobalIndex( row, localIdx );
+}
+
+/*
+template<>
+class MultidiagonalDeviceDependentCode< Devices::Host >
+{
+ public:
+
+ typedef Devices::Host Device;
+
+ template< typename Index >
+ __cuda_callable__
+ static Index getElementIndex( const Index rows,
+ const Index row,
+ const Index column )
+ {
+ return 2*row + column;
+ }
+
+ template< typename Vector,
+ typename Index,
+ typename ValuesType >
+ __cuda_callable__
+ static typename Vector::RealType rowVectorProduct( const Index rows,
+ const ValuesType& values,
+ const Index row,
+ const Vector& vector )
+ {
+ if( row == 0 )
+ return vector[ 0 ] * values[ 0 ] +
+ vector[ 1 ] * values[ 1 ];
+ Index i = 3 * row;
+ if( row == rows - 1 )
+ return vector[ row - 1 ] * values[ i - 1 ] +
+ vector[ row ] * values[ i ];
+ return vector[ row - 1 ] * values[ i - 1 ] +
+ vector[ row ] * values[ i ] +
+ vector[ row + 1 ] * values[ i + 1 ];
+ }
+
+ template< typename Real,
+ typename Index,
+ typename InVector,
+ typename OutVector >
+ static void vectorProduct( const Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >& matrix,
+ const InVector& inVector,
+ OutVector& outVector )
+ {
+#ifdef HAVE_OPENMP
+#pragma omp parallel for if( Devices::Host::isOMPEnabled() )
+#endif
+ for( Index row = 0; row < matrix.getRows(); row ++ )
+ outVector[ row ] = matrix.rowVectorProduct( row, inVector );
+ }
+};
+
+template<>
+class MultidiagonalDeviceDependentCode< Devices::Cuda >
+{
+ public:
+
+ typedef Devices::Cuda Device;
+
+ template< typename Index >
+ __cuda_callable__
+ static Index getElementIndex( const Index rows,
+ const Index row,
+ const Index column )
+ {
+ return ( column - row + 1 )*rows + row - 1;
+ }
+
+ template< typename Vector,
+ typename Index,
+ typename ValuesType >
+ __cuda_callable__
+ static typename Vector::RealType rowVectorProduct( const Index rows,
+ const ValuesType& values,
+ const Index row,
+ const Vector& vector )
+ {
+ if( row == 0 )
+ return vector[ 0 ] * values[ 0 ] +
+ vector[ 1 ] * values[ rows - 1 ];
+ Index i = row - 1;
+ if( row == rows - 1 )
+ return vector[ row - 1 ] * values[ i ] +
+ vector[ row ] * values[ i + rows ];
+ return vector[ row - 1 ] * values[ i ] +
+ vector[ row ] * values[ i + rows ] +
+ vector[ row + 1 ] * values[ i + 2*rows ];
+ }
+
+ template< typename Real,
+ typename Index,
+ typename InVector,
+ typename OutVector >
+ static void vectorProduct( const Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >& matrix,
+ const InVector& inVector,
+ OutVector& outVector )
+ {
+ MatrixVectorProductCuda( matrix, inVector, outVector );
+ }
+};
+ */
+
+} // namespace Matrices
+} // namespace TNL
diff --git a/src/TNL/Matrices/MultidiagonalMatrixRowView.h b/src/TNL/Matrices/MultidiagonalMatrixRowView.h
new file mode 100644
index 0000000000000000000000000000000000000000..68b5be55c1508994cf478ea9cb849e8c9710fcdc
--- /dev/null
+++ b/src/TNL/Matrices/MultidiagonalMatrixRowView.h
@@ -0,0 +1,59 @@
+/***************************************************************************
+ 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>
diff --git a/src/TNL/Matrices/MultidiagonalMatrixRowView.hpp b/src/TNL/Matrices/MultidiagonalMatrixRowView.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..349fbe8ea3eb430f41cf7481a8cfa27f91afdc36
--- /dev/null
+++ b/src/TNL/Matrices/MultidiagonalMatrixRowView.hpp
@@ -0,0 +1,75 @@
+/***************************************************************************
+ 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
diff --git a/src/TNL/Matrices/MultidiagonalMatrixView.h b/src/TNL/Matrices/MultidiagonalMatrixView.h
new file mode 100644
index 0000000000000000000000000000000000000000..addeb18b346340480c5c6eb41ec2b76893350b88
--- /dev/null
+++ b/src/TNL/Matrices/MultidiagonalMatrixView.h
@@ -0,0 +1,181 @@
+/***************************************************************************
+ 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>
diff --git a/src/TNL/Matrices/MultidiagonalMatrixView.hpp b/src/TNL/Matrices/MultidiagonalMatrixView.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..3d9b0237f908236da68b12f758b0f0fa51efb505
--- /dev/null
+++ b/src/TNL/Matrices/MultidiagonalMatrixView.hpp
@@ -0,0 +1,729 @@
+/***************************************************************************
+ MultidiagonalMatrixView.hpp - 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/Assert.h>
+#include <TNL/Matrices/MultidiagonalMatrixView.h>
+#include <TNL/Exceptions/NotImplementedError.h>
+
+namespace TNL {
+namespace Matrices {
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+MultidiagonalMatrixView()
+{
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+MultidiagonalMatrixView( const ValuesViewType& values,
+ const DiagonalsShiftsView& diagonalsShifts,
+ const IndexerType& indexer )
+: MatrixView< Real, Device, Index >( indexer.getRows(), indexer.getColumns(), values ),
+ diagonalsShifts( diagonalsShifts ),
+ indexer( indexer )
+{
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+auto
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+getView() -> ViewType
+{
+ return ViewType( this->values.getView(), indexer );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+auto
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+getConstView() const -> ConstViewType
+{
+ return ConstViewType( this->values.getConstView(), indexer );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+String
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+getSerializationType()
+{
+ return String( "Matrices::Multidiagonal< " ) +
+ TNL::getSerializationType< RealType >() + ", [any_device], " +
+ TNL::getSerializationType< IndexType >() + ", " +
+ ( RowMajorOrder ? "true" : "false" ) + ", [any_allocator] >";
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+String
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+getSerializationTypeVirtual() const
+{
+ return this->getSerializationType();
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+__cuda_callable__
+const Index&
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+getDiagonalsCount() const
+{
+ return this->diagonalsShifts.getSize();
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+ template< typename Vector >
+void
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+getCompressedRowLengths( Vector& rowLengths ) const
+{
+ rowLengths.setSize( this->getRows() );
+ rowLengths = 0;
+ auto rowLengths_view = rowLengths.getView();
+ auto fetch = [] __cuda_callable__ ( IndexType row, IndexType column, const RealType& value ) -> IndexType {
+ return ( value != 0.0 );
+ };
+ auto reduce = [] __cuda_callable__ ( IndexType& aux, const IndexType a ) {
+ aux += a;
+ };
+ auto keep = [=] __cuda_callable__ ( const IndexType rowIdx, const IndexType value ) mutable {
+ rowLengths_view[ rowIdx ] = value;
+ };
+ this->allRowsReduction( fetch, reduce, keep, 0 );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+Index
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+getNonemptyRowsCount() const
+{
+ return this->indexer.getNonemptyRowsCount();
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+Index
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+getRowLength( const IndexType row ) const
+{
+ return this->diagonalsShifts.getSize();
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+Index
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+getMaxRowLength() const
+{
+ return this->diagonalsShifts.getSize();
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+Index
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+getNumberOfNonzeroMatrixElements() const
+{
+ const auto values_view = this->values.getConstView();
+ auto fetch = [=] __cuda_callable__ ( const IndexType i ) -> IndexType {
+ return ( values_view[ i ] != 0.0 );
+ };
+ return Algorithms::Reduction< DeviceType >::reduce( this->values.getSize(), std::plus<>{}, fetch, 0 );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+ template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_ >
+bool
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+operator == ( const MultidiagonalMatrixView< Real_, Device_, Index_, RowMajorOrder_ >& matrix ) const
+{
+ if( RowMajorOrder == RowMajorOrder_ )
+ return this->values == matrix.values;
+ else
+ {
+ TNL_ASSERT( false, "TODO" );
+ }
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+ template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_ >
+bool
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+operator != ( const MultidiagonalMatrixView< Real_, Device_, Index_, RowMajorOrder_ >& matrix ) const
+{
+ return ! this->operator==( matrix );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+void
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+setValue( const RealType& v )
+{
+ this->values = v;
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+__cuda_callable__
+auto
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+getRow( const IndexType& rowIdx ) const -> const RowView
+{
+ return RowView( rowIdx, this->values.getView(), this->indexer );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+__cuda_callable__
+auto
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+getRow( const IndexType& rowIdx ) -> RowView
+{
+ return RowView( rowIdx, this->values.getView(), this->indexer );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+bool
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+setElement( const IndexType row, const IndexType column, const RealType& value )
+{
+ TNL_ASSERT_GE( row, 0, "" );
+ TNL_ASSERT_LT( row, this->getRows(), "" );
+ TNL_ASSERT_GE( column, 0, "" );
+ TNL_ASSERT_LT( column, this->getColumns(), "" );
+ if( abs( row - column ) > 1 )
+ {
+ std::stringstream msg;
+ msg << "Wrong matrix element coordinates ( " << row << ", " << column << " ) in tridiagonal matrix.";
+ throw std::logic_error( msg.str() );
+ }
+ this->values.setElement( this->getElementIndex( row, column ), value );
+ return true;
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+bool
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+addElement( const IndexType row,
+ const IndexType column,
+ const RealType& value,
+ const RealType& thisElementMultiplicator )
+{
+ TNL_ASSERT_GE( row, 0, "" );
+ TNL_ASSERT_LT( row, this->getRows(), "" );
+ TNL_ASSERT_GE( column, 0, "" );
+ TNL_ASSERT_LT( column, this->getColumns(), "" );
+ if( abs( row - column ) > 1 )
+ {
+ std::stringstream msg;
+ msg << "Wrong matrix element coordinates ( " << row << ", " << column << " ) in tridiagonal matrix.";
+ throw std::logic_error( msg.str() );
+ }
+ const Index i = this->getElementIndex( row, column );
+ this->values.setElement( i, thisElementMultiplicator * this->values.getElement( i ) + value );
+ return true;
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+Real
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+getElement( const IndexType row, const IndexType column ) const
+{
+ TNL_ASSERT_GE( row, 0, "" );
+ TNL_ASSERT_LT( row, this->getRows(), "" );
+ TNL_ASSERT_GE( column, 0, "" );
+ TNL_ASSERT_LT( column, this->getColumns(), "" );
+
+ if( abs( column - row ) > 1 )
+ return 0.0;
+ return this->values.getElement( this->getElementIndex( row, column ) );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >&
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+operator=( const MultidiagonalMatrixView& view )
+{
+ MatrixView< Real, Device, Index >::operator=( view );
+ this->diagonalsShifts.copy( view.diagonalsShifts );
+ this->hostDiagonalsShifts.copy( view.hostDiagonalsShifts );
+ this->indexer = view.indexer;
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+ template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >
+void
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+rowsReduction( IndexType first, IndexType last, Fetch& fetch, Reduce& reduce, Keep& keep, const FetchReal& zero_ ) const
+{
+ using Real_ = decltype( fetch( IndexType(), IndexType(), RealType() ) );
+ const auto values_view = this->values.getConstView();
+ const auto indexer = this->indexer;
+ const auto zero = zero_;
+ auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
+ Real_ sum( zero );
+ if( rowIdx == 0 )
+ {
+ reduce( sum, fetch( 0, 0, values_view[ indexer.getGlobalIndex( 0, 0 ) ] ) );
+ reduce( sum, fetch( 0, 1, values_view[ indexer.getGlobalIndex( 0, 1 ) ] ) );
+ keep( 0, sum );
+ return;
+ }
+ if( rowIdx + 1 < indexer.getColumns() )
+ {
+ reduce( sum, fetch( rowIdx, rowIdx - 1, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ] ) );
+ reduce( sum, fetch( rowIdx, rowIdx, values_view[ indexer.getGlobalIndex( rowIdx, 1 ) ] ) );
+ reduce( sum, fetch( rowIdx, rowIdx + 1, values_view[ indexer.getGlobalIndex( rowIdx, 2 ) ] ) );
+ keep( rowIdx, sum );
+ return;
+ }
+ if( rowIdx < indexer.getColumns() )
+ {
+ reduce( sum, fetch( rowIdx, rowIdx - 1, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ] ) );
+ reduce( sum, fetch( rowIdx, rowIdx, values_view[ indexer.getGlobalIndex( rowIdx, 1 ) ] ) );
+ keep( rowIdx, sum );
+ }
+ else
+ {
+ keep( rowIdx, fetch( rowIdx, rowIdx - 1, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ] ) );
+ }
+ };
+ Algorithms::ParallelFor< DeviceType >::exec( first, last, f );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+ template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >
+void
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+allRowsReduction( Fetch& fetch, Reduce& reduce, Keep& keep, const FetchReal& zero ) const
+{
+ this->rowsReduction( 0, this->indexer.getNonEmptyRowsCount(), fetch, reduce, keep, zero );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+ template< typename Function >
+void
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+forRows( IndexType first, IndexType last, Function& function ) const
+{
+ const auto values_view = this->values.getConstView();
+ const auto indexer = this->indexer;
+ auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
+ if( rowIdx == 0 )
+ {
+ function( 0, 0, 0, values_view[ indexer.getGlobalIndex( 0, 0 ) ] );
+ function( 0, 1, 1, values_view[ indexer.getGlobalIndex( 0, 1 ) ] );
+ }
+ else if( rowIdx + 1 < indexer.getColumns() )
+ {
+ function( rowIdx, 0, rowIdx - 1, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ] );
+ function( rowIdx, 1, rowIdx, values_view[ indexer.getGlobalIndex( rowIdx, 1 ) ] );
+ function( rowIdx, 2, rowIdx + 1, values_view[ indexer.getGlobalIndex( rowIdx, 2 ) ] );
+ }
+ else if( rowIdx < indexer.getColumns() )
+ {
+ function( rowIdx, 0, rowIdx - 1, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ] );
+ function( rowIdx, 1, rowIdx, values_view[ indexer.getGlobalIndex( rowIdx, 1 ) ] );
+ }
+ else
+ function( rowIdx, 0, rowIdx, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ] );
+ };
+ Algorithms::ParallelFor< DeviceType >::exec( first, last, f );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+ template< typename Function >
+void
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+forRows( IndexType first, IndexType last, Function& function )
+{
+ auto values_view = this->values.getView();
+ const auto indexer = this->indexer;
+ auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
+ if( rowIdx == 0 )
+ {
+ function( 0, 0, 0, values_view[ indexer.getGlobalIndex( 0, 0 ) ] );
+ function( 0, 1, 1, values_view[ indexer.getGlobalIndex( 0, 1 ) ] );
+ }
+ else if( rowIdx + 1 < indexer.getColumns() )
+ {
+ function( rowIdx, 0, rowIdx - 1, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ] );
+ function( rowIdx, 1, rowIdx, values_view[ indexer.getGlobalIndex( rowIdx, 1 ) ] );
+ function( rowIdx, 2, rowIdx + 1, values_view[ indexer.getGlobalIndex( rowIdx, 2 ) ] );
+ }
+ else if( rowIdx < indexer.getColumns() )
+ {
+ function( rowIdx, 0, rowIdx - 1, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ] );
+ function( rowIdx, 1, rowIdx, values_view[ indexer.getGlobalIndex( rowIdx, 1 ) ] );
+ }
+ else
+ function( rowIdx, 0, rowIdx, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ] );
+ };
+ Algorithms::ParallelFor< DeviceType >::exec( first, last, f );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+ template< typename Function >
+void
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+forAllRows( Function& function ) const
+{
+ this->forRows( 0, this->indxer.getNonEmptyRowsCount(), function );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+ template< typename Function >
+void
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+forAllRows( Function& function )
+{
+ this->forRows( 0, this->indexer.getNonEmptyRowsCount(), function );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+template< typename Vector >
+__cuda_callable__
+typename Vector::RealType
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+rowVectorProduct( const IndexType row, const Vector& vector ) const
+{
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+ template< typename InVector,
+ typename OutVector >
+void
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+vectorProduct( const InVector& inVector, OutVector& outVector ) const
+{
+ TNL_ASSERT_EQ( this->getColumns(), inVector.getSize(), "Matrix columns do not fit with input vector." );
+ TNL_ASSERT_EQ( this->getRows(), outVector.getSize(), "Matrix rows do not fit with output vector." );
+
+ const auto inVectorView = inVector.getConstView();
+ auto outVectorView = outVector.getView();
+ const auto valuesView = this->values.getConstView();
+ auto fetch = [=] __cuda_callable__ ( const IndexType& row, const IndexType& column, const RealType& value ) -> RealType {
+ return value * inVectorView[ column ];
+ };
+ auto reduction = [] __cuda_callable__ ( RealType& sum, const RealType& value ) {
+ sum += value;
+ };
+ auto keeper = [=] __cuda_callable__ ( IndexType row, const RealType& value ) mutable {
+ outVectorView[ row ] = value;
+ };
+ this->allRowsReduction( fetch, reduction, keeper, ( RealType ) 0.0 );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+ template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_ >
+void
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+addMatrix( const MultidiagonalMatrixView< Real_, Device_, Index_, RowMajorOrder_ >& matrix,
+ const RealType& matrixMultiplicator,
+ const RealType& thisMatrixMultiplicator )
+{
+ TNL_ASSERT_EQ( this->getRows(), matrix.getRows(), "Matrices rows are not equal." );
+ TNL_ASSERT_EQ( this->getColumns(), matrix.getColumns(), "Matrices columns are not equal." );
+
+ if( RowMajorOrder == RowMajorOrder_ )
+ {
+ if( thisMatrixMultiplicator == 1.0 )
+ this->values += matrixMultiplicator * matrix.getValues();
+ else
+ this->values = thisMatrixMultiplicator * this->values + matrixMultiplicator * matrix.getValues();
+ }
+ else
+ {
+ const auto matrix_view = matrix;
+ const auto matrixMult = matrixMultiplicator;
+ const auto thisMult = thisMatrixMultiplicator;
+ auto add0 = [=] __cuda_callable__ ( const IndexType& rowIdx, const IndexType& localIdx, const IndexType& column, Real& value ) mutable {
+ value = matrixMult * matrix.getValues()[ matrix.getIndexer().getGlobalIndex( rowIdx, localIdx ) ];
+ };
+ auto add1 = [=] __cuda_callable__ ( const IndexType& rowIdx, const IndexType& localIdx, const IndexType& column, Real& value ) mutable {
+ value += matrixMult * matrix.getValues()[ matrix.getIndexer().getGlobalIndex( rowIdx, localIdx ) ];
+ };
+ auto addGen = [=] __cuda_callable__ ( const IndexType& rowIdx, const IndexType& localIdx, const IndexType& column, Real& value ) mutable {
+ value = thisMult * value + matrixMult * matrix.getValues()[ matrix.getIndexer().getGlobalIndex( rowIdx, localIdx ) ];
+ };
+ if( thisMult == 0.0 )
+ this->forAllRows( add0 );
+ else if( thisMult == 1.0 )
+ this->forAllRows( add1 );
+ else
+ this->forAllRows( addGen );
+ }
+}
+
+#ifdef HAVE_CUDA
+/*template< typename Real,
+ typename Real2,
+ typename Index,
+ typename Index2 >
+__global__ void MultidiagonalTranspositionCudaKernel( const Multidiagonal< Real2, Devices::Cuda, Index2 >* inMatrix,
+ Multidiagonal< Real, Devices::Cuda, Index >* outMatrix,
+ const Real matrixMultiplicator,
+ const Index gridIdx )
+{
+ const Index rowIdx = ( gridIdx * Cuda::getMaxGridSize() + blockIdx.x ) * blockDim.x + threadIdx.x;
+ if( rowIdx < inMatrix->getRows() )
+ {
+ if( rowIdx > 0 )
+ outMatrix->setElementFast( rowIdx-1,
+ rowIdx,
+ matrixMultiplicator * inMatrix->getElementFast( rowIdx, rowIdx-1 ) );
+ outMatrix->setElementFast( rowIdx,
+ rowIdx,
+ matrixMultiplicator * inMatrix->getElementFast( rowIdx, rowIdx ) );
+ if( rowIdx < inMatrix->getRows()-1 )
+ outMatrix->setElementFast( rowIdx+1,
+ rowIdx,
+ matrixMultiplicator * inMatrix->getElementFast( rowIdx, rowIdx+1 ) );
+ }
+}*/
+#endif
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+ template< typename Real2, typename Index2 >
+void
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+getTransposition( const MultidiagonalMatrixView< Real2, Device, Index2 >& matrix,
+ const RealType& matrixMultiplicator )
+{
+ TNL_ASSERT( this->getRows() == matrix.getRows(),
+ std::cerr << "This matrix rows: " << this->getRows() << std::endl
+ << "That matrix rows: " << matrix.getRows() << std::endl );
+ if( std::is_same< Device, Devices::Host >::value )
+ {
+ const IndexType& rows = matrix.getRows();
+ for( IndexType i = 1; i < rows; i++ )
+ {
+ RealType aux = matrix. getElement( i, i - 1 );
+ this->setElement( i, i - 1, matrix.getElement( i - 1, i ) );
+ this->setElement( i, i, matrix.getElement( i, i ) );
+ this->setElement( i - 1, i, aux );
+ }
+ }
+ if( std::is_same< Device, Devices::Cuda >::value )
+ {
+#ifdef HAVE_CUDA
+ /*Multidiagonal* kernel_this = Cuda::passToDevice( *this );
+ typedef Multidiagonal< Real2, Device, Index2 > InMatrixType;
+ InMatrixType* kernel_inMatrix = Cuda::passToDevice( matrix );
+ dim3 cudaBlockSize( 256 ), cudaGridSize( Cuda::getMaxGridSize() );
+ const IndexType cudaBlocks = roundUpDivision( matrix.getRows(), cudaBlockSize.x );
+ const IndexType cudaGrids = roundUpDivision( cudaBlocks, Cuda::getMaxGridSize() );
+ for( IndexType gridIdx = 0; gridIdx < cudaGrids; gridIdx++ )
+ {
+ if( gridIdx == cudaGrids - 1 )
+ cudaGridSize.x = cudaBlocks % Cuda::getMaxGridSize();
+ MultidiagonalTranspositionCudaKernel<<< cudaGridSize, cudaBlockSize >>>
+ ( kernel_inMatrix,
+ kernel_this,
+ matrixMultiplicator,
+ gridIdx );
+ }
+ Cuda::freeFromDevice( kernel_this );
+ Cuda::freeFromDevice( kernel_inMatrix );
+ TNL_CHECK_CUDA_DEVICE;*/
+#endif
+ }
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+ template< typename Vector1, typename Vector2 >
+__cuda_callable__
+void
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+performSORIteration( const Vector1& b,
+ const IndexType row,
+ Vector2& x,
+ const RealType& omega ) const
+{
+ RealType sum( 0.0 );
+ if( row > 0 )
+ sum += this->getElementFast( row, row - 1 ) * x[ row - 1 ];
+ if( row < this->getColumns() - 1 )
+ sum += this->getElementFast( row, row + 1 ) * x[ row + 1 ];
+ x[ row ] = ( 1.0 - omega ) * x[ row ] + omega / this->getElementFast( row, row ) * ( b[ row ] - sum );
+}
+
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+void MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::save( File& file ) const
+{
+ MatrixView< Real, Device, Index >::save( file );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+void
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+save( const String& fileName ) const
+{
+ Object::save( fileName );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+void MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::print( std::ostream& str ) const
+{
+ for( IndexType row = 0; row < this->getRows(); row++ )
+ {
+ str <<"Row: " << row << " -> ";
+ for( IndexType column = row - 1; column < row + 2; column++ )
+ if( column >= 0 && column < this->columns )
+ str << " Col:" << column << "->" << this->getElement( row, column ) << "\t";
+ str << std::endl;
+ }
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+__cuda_callable__
+auto
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+getIndexer() const -> const IndexerType&
+{
+ return this->indexer;
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+__cuda_callable__
+auto
+MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+getIndexer() -> IndexerType&
+{
+ return this->indexer;
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ bool RowMajorOrder >
+__cuda_callable__
+Index MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
+getElementIndex( const IndexType row, const IndexType column ) const
+{
+ IndexType localIdx = column - row;
+ if( row > 0 )
+ localIdx++;
+
+ TNL_ASSERT_GE( localIdx, 0, "" );
+ TNL_ASSERT_LT( localIdx, 3, "" );
+
+ return this->indexer.getGlobalIndex( row, localIdx );
+}
+
+} // namespace Matrices
+} // namespace TNL
diff --git a/src/TNL/Matrices/TridiagonalMatrixView.h b/src/TNL/Matrices/TridiagonalMatrixView.h
index 29006279390150cbd80a469c8ffed124ca67be17..128b4849476917bd0f7417a5f881b3d72c6c7be7 100644
--- a/src/TNL/Matrices/TridiagonalMatrixView.h
+++ b/src/TNL/Matrices/TridiagonalMatrixView.h
@@ -59,9 +59,6 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
virtual String getSerializationTypeVirtual() const;
- void setDimensions( const IndexType rows,
- const IndexType columns );
-
template< typename Vector >
void getCompressedRowLengths( Vector& rowLengths ) const;
diff --git a/src/TNL/Matrices/details/MultidiagonalMatrixIndexer.h b/src/TNL/Matrices/details/MultidiagonalMatrixIndexer.h
new file mode 100644
index 0000000000000000000000000000000000000000..0f0436d748c9713e106964398387bd49ede6042e
--- /dev/null
+++ b/src/TNL/Matrices/details/MultidiagonalMatrixIndexer.h
@@ -0,0 +1,106 @@
+/***************************************************************************
+ MultidiagonalMatrixIndexer.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
+
+namespace TNL {
+ namespace Matrices {
+ namespace details {
+
+template< typename Index,
+ bool RowMajorOrder >
+class MultidiagonalMatrixIndexer
+{
+ public:
+
+ using IndexType = Index;
+
+ static constexpr bool getRowMajorOrder() { return RowMajorOrder; };
+
+ __cuda_callable__
+ MultidiagonalMatrixIndexer()
+ : rows( 0 ), columns( 0 ), nonemptyRows( 0 ){};
+
+ __cuda_callable__
+ MultidiagonalMatrixIndexer( const IndexType& rows,
+ const IndexType& columns,
+ const IndexType& diagonals,
+ const IndexType& nonemptyRows )
+ : rows( rows ),
+ columns( columns ),
+ diagonals( diagonals ),
+ nonemptyRows( nonemptyRows ) {};
+
+ __cuda_callable__
+ MultidiagonalMatrixIndexer( const MultidiagonalMatrixIndexer& indexer )
+ : rows( indexer.rows ),
+ columns( indexer.columns ),
+ diagonals( indexer.diagonals ),
+ nonemptyRows( indexer.nonemptyRows ) {};
+
+ void set( const IndexType& rows,
+ const IndexType& columns,
+ const IndexType& diagonals,
+ const IndexType& nonemptyRows )
+ {
+ this->rows = rows;
+ this->columns = columns;
+ this->diagonals = diagonals;
+ this->nonemptyRows = nonemptyRows;
+ };
+
+ /*__cuda_callable__
+ IndexType getRowSize( const IndexType rowIdx ) const
+ {
+ if( rowIdx == 0 )
+ return 2;
+ if( columns <= rows )
+ {
+ if( rowIdx == columns - 1 )
+ return 2;
+ if( rowIdx == columns )
+ return 1;
+ }
+ return 3;
+ };*/
+
+ __cuda_callable__
+ const IndexType& getRows() const { return this->rows; };
+
+ __cuda_callable__
+ const IndexType& getColumns() const { return this->columns; };
+
+ __cuda_callable__
+ const IndexType& getNonemptyRowsCount() const { return this->nonemptyRows; };
+
+ __cuda_callable__
+ IndexType getStorageSize() const { return diagonals * this->nonemptyRows; };
+
+ __cuda_callable__
+ IndexType getGlobalIndex( const Index rowIdx, const Index localIdx ) const
+ {
+ TNL_ASSERT_GE( localIdx, 0, "" );
+ TNL_ASSERT_LT( localIdx, diagonals, "" );
+ TNL_ASSERT_GE( rowIdx, 0, "" );
+ TNL_ASSERT_LT( rowIdx, this->rows, "" );
+
+ if( RowMajorOrder )
+ return diagonals * rowIdx + localIdx;
+ else
+ return localIdx * nonemptyRows + rowIdx;
+ };
+
+ protected:
+
+ IndexType rows, columns, diagonals, nonemptyRows;
+};
+ } //namespace details
+ } // namespace Materices
+} // namespace TNL
diff --git a/src/UnitTests/Matrices/CMakeLists.txt b/src/UnitTests/Matrices/CMakeLists.txt
index 4b95380c4477b33184e6c80946d9f0aebd95e04c..28749540556121c826e5de3457ddad9f3996201f 100644
--- a/src/UnitTests/Matrices/CMakeLists.txt
+++ b/src/UnitTests/Matrices/CMakeLists.txt
@@ -13,8 +13,8 @@ IF( BUILD_CUDA )
CUDA_ADD_EXECUTABLE( TridiagonalMatrixTest TridiagonalMatrixTest.cu OPTIONS ${CXX_TESTS_FLAGS} )
TARGET_LINK_LIBRARIES( TridiagonalMatrixTest ${GTEST_BOTH_LIBRARIES} )
-# CUDA_ADD_EXECUTABLE( MultidiagonalMatrixTest MultidiagonalMatrixTest.cu OPTIONS ${CXX_TESTS_FLAGS} )
-# TARGET_LINK_LIBRARIES( MultidiagonalMatrixTest ${GTEST_BOTH_LIBRARIES} )
+ CUDA_ADD_EXECUTABLE( MultidiagonalMatrixTest MultidiagonalMatrixTest.cu OPTIONS ${CXX_TESTS_FLAGS} )
+ TARGET_LINK_LIBRARIES( MultidiagonalMatrixTest ${GTEST_BOTH_LIBRARIES} )
CUDA_ADD_EXECUTABLE( SparseMatrixTest_CSR_segments SparseMatrixTest_CSR_segments.cu OPTIONS ${CXX_TESTS_FLAGS} )
TARGET_LINK_LIBRARIES( SparseMatrixTest_CSR_segments ${GTEST_BOTH_LIBRARIES} )
@@ -42,9 +42,9 @@ ELSE( BUILD_CUDA )
TARGET_COMPILE_OPTIONS( TridiagonalMatrixTest PRIVATE ${CXX_TESTS_FLAGS} )
TARGET_LINK_LIBRARIES( TridiagonalMatrixTest ${GTEST_BOTH_LIBRARIES} )
-# ADD_EXECUTABLE( MultidiagonalMatrixTest MultidiagonalMatrixTest.cpp )
-# TARGET_COMPILE_OPTIONS( MultidiagonalMatrixTest PRIVATE ${CXX_TESTS_FLAGS} )
-# TARGET_LINK_LIBRARIES( MultidiagonalMatrixTest ${GTEST_BOTH_LIBRARIES} )
+ ADD_EXECUTABLE( MultidiagonalMatrixTest MultidiagonalMatrixTest.cpp )
+ TARGET_COMPILE_OPTIONS( MultidiagonalMatrixTest PRIVATE ${CXX_TESTS_FLAGS} )
+ TARGET_LINK_LIBRARIES( MultidiagonalMatrixTest ${GTEST_BOTH_LIBRARIES} )
ADD_EXECUTABLE( SparseMatrixTest_CSR_segments SparseMatrixTest_CSR_segments.cpp )
TARGET_COMPILE_OPTIONS( SparseMatrixTest_CSR_segments PRIVATE ${CXX_TESTS_FLAGS} )
@@ -65,7 +65,7 @@ ADD_TEST( SparseMatrixCopyTest ${EXECUTABLE_OUTPUT_PATH}/SparseMatrixCopyTest${C
ADD_TEST( SparseMatrixTest ${EXECUTABLE_OUTPUT_PATH}/SparseMatrixTest${CMAKE_EXECUTABLE_SUFFIX} )
ADD_TEST( DenseMatrixTest ${EXECUTABLE_OUTPUT_PATH}/DenseMatrixTest${CMAKE_EXECUTABLE_SUFFIX} )
ADD_TEST( TridiagonalMatrixTest ${EXECUTABLE_OUTPUT_PATH}/TridiagonalMatrixTest${CMAKE_EXECUTABLE_SUFFIX} )
-#ADD_TEST( MultidiagonalMatrixTest ${EXECUTABLE_OUTPUT_PATH}/MultidiagonalMatrixTest${CMAKE_EXECUTABLE_SUFFIX} )
+ADD_TEST( MultidiagonalMatrixTest ${EXECUTABLE_OUTPUT_PATH}/MultidiagonalMatrixTest${CMAKE_EXECUTABLE_SUFFIX} )
ADD_TEST( SparseMatrixTest_CSR_segments ${EXECUTABLE_OUTPUT_PATH}/SparseMatrixTest_CSR_segments${CMAKE_EXECUTABLE_SUFFIX} )
ADD_TEST( SparseMatrixTest_Ellpack_segments ${EXECUTABLE_OUTPUT_PATH}/SparseMatrixTest_Ellpack_segments${CMAKE_EXECUTABLE_SUFFIX} )
diff --git a/src/UnitTests/Matrices/MultidiagonalMatrixTest.cpp b/src/UnitTests/Matrices/MultidiagonalMatrixTest.cpp
index 73406d0dfc49879c507fded1ff3a06b7cad78639..639f1964086784bfdb174443a5f2554b703e511b 100644
--- a/src/UnitTests/Matrices/MultidiagonalMatrixTest.cpp
+++ b/src/UnitTests/Matrices/MultidiagonalMatrixTest.cpp
@@ -1,7 +1,7 @@
/***************************************************************************
MultidiagonalMatrixTest.cpp - description
-------------------
- begin : Jan 9, 2020
+ begin : Jan 8, 2020
copyright : (C) 2020 by Tomas Oberhuber et al.
email : tomas.oberhuber@fjfi.cvut.cz
***************************************************************************/
diff --git a/src/UnitTests/Matrices/MultidiagonalMatrixTest.cu b/src/UnitTests/Matrices/MultidiagonalMatrixTest.cu
index e3dab545cb5bf9494e3f9197893275b6b88c7d49..53541edbd003d084e1b50b742472beec086d87fb 100644
--- a/src/UnitTests/Matrices/MultidiagonalMatrixTest.cu
+++ b/src/UnitTests/Matrices/MultidiagonalMatrixTest.cu
@@ -1,7 +1,7 @@
/***************************************************************************
MultidiagonalMatrixTest.cu - description
-------------------
- begin : Jan 9, 2020
+ begin : Jan 8, 2020
copyright : (C) 2020 by Tomas Oberhuber et al.
email : tomas.oberhuber@fjfi.cvut.cz
***************************************************************************/
diff --git a/src/UnitTests/Matrices/MultidiagonalMatrixTest.h b/src/UnitTests/Matrices/MultidiagonalMatrixTest.h
index abe6b64c555fcb54068a5e53df968cfe8e11d10a..cb9916e4c939d11483c9229ef03da6e5eb661483 100644
--- a/src/UnitTests/Matrices/MultidiagonalMatrixTest.h
+++ b/src/UnitTests/Matrices/MultidiagonalMatrixTest.h
@@ -8,6 +8,7 @@
/* See Copyright Notice in tnl/Copyright */
+#include <sstream>
#include <TNL/Devices/Host.h>
#include <TNL/Matrices/Matrix.h>
#include <TNL/Matrices/Multidiagonal.h>
@@ -33,55 +34,110 @@ static const char* TEST_FILE_NAME = "test_MultidiagonalMatrixTest.tnl";
void test_GetSerializationType()
{
- EXPECT_EQ( ( TNL::Matrices::Multidiagonal< float, TNL::Devices::Host, int, true >::getSerializationType() ), TNL::String( "Matrices::Multidiagonal< float, [any_device], int, true, [any_allocator] >" ) );
- EXPECT_EQ( ( TNL::Matrices::Multidiagonal< int, TNL::Devices::Host, int, true >::getSerializationType() ), TNL::String( "Matrices::Multidiagonal< int, [any_device], int, true, [any_allocator] >" ) );
- EXPECT_EQ( ( TNL::Matrices::Multidiagonal< float, TNL::Devices::Cuda, int, true >::getSerializationType() ), TNL::String( "Matrices::Multidiagonal< float, [any_device], int, true, [any_allocator] >" ) );
- EXPECT_EQ( ( TNL::Matrices::Multidiagonal< int, TNL::Devices::Cuda, int, true >::getSerializationType() ), TNL::String( "Matrices::Multidiagonal< int, [any_device], int, true, [any_allocator] >" ) );
- EXPECT_EQ( ( TNL::Matrices::Multidiagonal< float, TNL::Devices::Host, int, false >::getSerializationType() ), TNL::String( "Matrices::Multidiagonal< float, [any_device], int, false, [any_allocator] >" ) );
- EXPECT_EQ( ( TNL::Matrices::Multidiagonal< int, TNL::Devices::Host, int, false >::getSerializationType() ), TNL::String( "Matrices::Multidiagonal< int, [any_device], int, false, [any_allocator] >" ) );
- EXPECT_EQ( ( TNL::Matrices::Multidiagonal< float, TNL::Devices::Cuda, int, false >::getSerializationType() ), TNL::String( "Matrices::Multidiagonal< float, [any_device], int, false, [any_allocator] >" ) );
- EXPECT_EQ( ( TNL::Matrices::Multidiagonal< int, TNL::Devices::Cuda, int, false >::getSerializationType() ), TNL::String( "Matrices::Multidiagonal< int, [any_device], int, false, [any_allocator] >" ) );
+ EXPECT_EQ( ( TNL::Matrices::Multidiagonal< float, TNL::Devices::Host, int, true >::getSerializationType() ), TNL::String( "Matrices::Multidiagonal< float, [any_device], int, true, [any_allocator], [any_allocator] >" ) );
+ EXPECT_EQ( ( TNL::Matrices::Multidiagonal< int, TNL::Devices::Host, int, true >::getSerializationType() ), TNL::String( "Matrices::Multidiagonal< int, [any_device], int, true, [any_allocator], [any_allocator] >" ) );
+ EXPECT_EQ( ( TNL::Matrices::Multidiagonal< float, TNL::Devices::Cuda, int, true >::getSerializationType() ), TNL::String( "Matrices::Multidiagonal< float, [any_device], int, true, [any_allocator], [any_allocator] >" ) );
+ EXPECT_EQ( ( TNL::Matrices::Multidiagonal< int, TNL::Devices::Cuda, int, true >::getSerializationType() ), TNL::String( "Matrices::Multidiagonal< int, [any_device], int, true, [any_allocator], [any_allocator] >" ) );
+ EXPECT_EQ( ( TNL::Matrices::Multidiagonal< float, TNL::Devices::Host, int, false >::getSerializationType() ), TNL::String( "Matrices::Multidiagonal< float, [any_device], int, false, [any_allocator], [any_allocator] >" ) );
+ EXPECT_EQ( ( TNL::Matrices::Multidiagonal< int, TNL::Devices::Host, int, false >::getSerializationType() ), TNL::String( "Matrices::Multidiagonal< int, [any_device], int, false, [any_allocator], [any_allocator] >" ) );
+ EXPECT_EQ( ( TNL::Matrices::Multidiagonal< float, TNL::Devices::Cuda, int, false >::getSerializationType() ), TNL::String( "Matrices::Multidiagonal< float, [any_device], int, false, [any_allocator], [any_allocator] >" ) );
+ EXPECT_EQ( ( TNL::Matrices::Multidiagonal< int, TNL::Devices::Cuda, int, false >::getSerializationType() ), TNL::String( "Matrices::Multidiagonal< int, [any_device], int, false, [any_allocator], [any_allocator] >" ) );
}
template< typename Matrix >
void test_SetDimensions()
{
- using RealType = typename Matrix::RealType;
- using DeviceType = typename Matrix::DeviceType;
- using IndexType = typename Matrix::IndexType;
+ using RealType = typename Matrix::RealType;
+ using DeviceType = typename Matrix::DeviceType;
+ using IndexType = typename Matrix::IndexType;
+ using DiagonalsShiftsType = typename Matrix::DiagonalsShiftsType;
- const IndexType rows = 9;
- const IndexType cols = 8;
+ const IndexType rows = 9;
+ const IndexType cols = 8;
+ const DiagonalsShiftsType diagonalsShifts{ -3, -1, 0, 2, 4 };
- Matrix m;
- m.setDimensions( rows, cols );
+ Matrix m;
+ m.setDimensions( rows, cols, diagonalsShifts );
- EXPECT_EQ( m.getRows(), 9 );
- EXPECT_EQ( m.getColumns(), 8 );
+ EXPECT_EQ( m.getRows(), 9 );
+ EXPECT_EQ( m.getColumns(), 8 );
}
+
template< typename Matrix1, typename Matrix2 >
void test_SetLike()
{
- using RealType = typename Matrix1::RealType;
- using DeviceType = typename Matrix1::DeviceType;
- using IndexType = typename Matrix1::IndexType;
+ using RealType = typename Matrix1::RealType;
+ using DeviceType = typename Matrix1::DeviceType;
+ using IndexType = typename Matrix1::IndexType;
+ using DiagonalsShiftsType = typename Matrix1::DiagonalsShiftsType;
- const IndexType rows = 8;
- const IndexType cols = 7;
+ const IndexType rows = 8;
+ const IndexType cols = 7;
+ const DiagonalsShiftsType diagonalsShifts{ -3, -1, 0, 2, 4 };
- Matrix1 m1;
- m1.reset();
- m1.setDimensions( rows + 1, cols + 2 );
+ Matrix1 m1;
+ m1.setDimensions( rows + 1, cols + 2, diagonalsShifts );
- Matrix2 m2;
- m2.reset();
- m2.setDimensions( rows, cols );
+ Matrix2 m2;
+ m2.setDimensions( rows, cols, diagonalsShifts );
+
+ m1.setLike( m2 );
+
+ EXPECT_EQ( m1.getRows(), m2.getRows() );
+ EXPECT_EQ( m1.getColumns(), m2.getColumns() );
+}
+
+template< typename Matrix >
+void test_GetNonemptyRowsCount()
+{
+ using RealType = typename Matrix::RealType;
+ using DeviceType = typename Matrix::DeviceType;
+ using IndexType = typename Matrix::IndexType;
+ using DiagonalsShiftsType = typename Matrix::DiagonalsShiftsType;
+
+ /*
+ * Sets up the following 5x8 matrix:
+ *
+ * / 1 0 0 1 0 1 0 0 \
+ * | 0 1 0 0 1 0 1 0 |
+ * | 1 0 1 0 0 1 0 1 |
+ * | 0 1 0 1 0 0 1 0 |
+ * \ 0 0 1 0 1 0 0 1 /
+ */
+ Matrix m1( 5, 8, DiagonalsShiftsType({ -2, 0, 3, 5 }) );
+ m1.setValue( 1.0 );
+ EXPECT_EQ( m1.getNonemptyRowsCount(), 5 );
+
+ /*
+ * Sets up the following 5x5 matrix:
+ *
+ * / 1 0 0 1 0 \
+ * | 0 1 0 0 1 |
+ * | 1 0 1 0 0 |
+ * | 0 1 0 1 0 |
+ * \ 0 0 1 0 1 /
+ */
+ Matrix m2( 5, 5, DiagonalsShiftsType({ -2, 0, 3, 5 }) );
+ m2.setValue( 1.0 );
+ EXPECT_EQ( m2.getNonemptyRowsCount(), 5 );
- m1.setLike( m2 );
+ /*
+ * Sets up the following 8x5 matrix:
+ *
+ * / 1 0 0 1 0 \
+ * | 0 1 0 0 1 |
+ * | 1 0 1 0 0 |
+ * | 0 1 0 1 0 |
+ * | 0 0 1 0 1 |
+ * | 0 0 0 1 0 |
+ * | 0 0 0 0 1 |
+ * \ 0 0 0 0 0 /
+ */
+ Matrix m3( 8, 5, DiagonalsShiftsType({ -2, 0, 3, 5 }) );
+ m3.setValue( 1.0 );
+ EXPECT_EQ( m3.getNonemptyRowsCount(), 7 );
- EXPECT_EQ( m1.getRows(), m2.getRows() );
- EXPECT_EQ( m1.getColumns(), m2.getColumns() );
}
template< typename Matrix >
@@ -90,463 +146,470 @@ void test_GetCompressedRowLengths()
using RealType = typename Matrix::RealType;
using DeviceType = typename Matrix::DeviceType;
using IndexType = typename Matrix::IndexType;
+ using DiagonalsShiftsType = typename Matrix::DiagonalsShiftsType;
const IndexType rows = 10;
const IndexType cols = 11;
- Matrix m( rows, cols );
+ Matrix m( rows, cols );
- // Insert values into the rows.
- RealType value = 1;
+ // Insert values into the rows.
+ RealType value = 1;
- for( IndexType i = 0; i < 3; i++ ) // 0th row
- m.setElement( 0, i, value++ );
+ for( IndexType i = 0; i < 2; i++ ) // 0th row -> 2 elements
+ m.setElement( 0, i, value++ );
- for( IndexType i = 0; i < 3; i++ ) // 1st row
- m.setElement( 1, i, value++ );
+ for( IndexType i = 0; i < 3; i++ ) // 1st row -> 3 elements
+ m.setElement( 1, i, value++ );
- for( IndexType i = 0; i < 1; i++ ) // 2nd row
- m.setElement( 2, i, value++ );
+ for( IndexType i = 1; i < 3; i++ ) // 2nd row -> 2 elements
+ m.setElement( 2, i, value++ );
- for( IndexType i = 0; i < 2; i++ ) // 3rd row
- m.setElement( 3, i, value++ );
+ for( IndexType i = 2; i < 5; i++ ) // 3rd row -> 3 elements
+ m.setElement( 3, i, value++ );
- for( IndexType i = 0; i < 3; i++ ) // 4th row
- m.setElement( 4, i, value++ );
+ for( IndexType i = 3; i < 6; i++ ) // 4th row -> 3 elements
+ m.setElement( 4, i, value++ );
- for( IndexType i = 0; i < 4; i++ ) // 5th row
- m.setElement( 5, i, value++ );
+ for( IndexType i = 4; i < 6; i++ ) // 5th row -> 2 elements
+ m.setElement( 5, i, value++ );
- for( IndexType i = 0; i < 5; i++ ) // 6th row
- m.setElement( 6, i, value++ );
+ for( IndexType i = 5; i < 8; i++ ) // 6th row -> 3 elements
+ m.setElement( 6, i, value++ );
- for( IndexType i = 0; i < 6; i++ ) // 7th row
- m.setElement( 7, i, value++ );
+ for( IndexType i = 6; i < 8; i++ ) // 7th row -> 2 elements
+ m.setElement( 7, i, value++ );
- for( IndexType i = 0; i < 7; i++ ) // 8th row
- m.setElement( 8, i, value++ );
+ for( IndexType i = 7; i < 10; i++ ) // 8th row -> 3 elements
+ m.setElement( 8, i, value++ );
- for( IndexType i = 0; i < 8; i++ ) // 9th row
- m.setElement( 9, i, value++ );
+ for( IndexType i = 8; i < 11; i++ ) // 9th row -> 3 elements
+ m.setElement( 9, i, value++ );
- typename Matrix::CompressedRowLengthsVector rowLengths;
+ typename Matrix::CompressedRowLengthsVector rowLengths( rows );
rowLengths = 0;
m.getCompressedRowLengths( rowLengths );
- typename Matrix::CompressedRowLengthsVector correctRowLengths{ 3, 3, 1, 2, 3, 4, 5, 6, 7, 8 };
+ typename Matrix::CompressedRowLengthsVector correctRowLengths{ 2, 3, 2, 3, 3, 2, 3, 2, 3, 3 };
EXPECT_EQ( rowLengths, correctRowLengths );
}
template< typename Matrix >
void test_GetRowLength()
{
- using RealType = typename Matrix::RealType;
- using DeviceType = typename Matrix::DeviceType;
- using IndexType = typename Matrix::IndexType;
+ using RealType = typename Matrix::RealType;
+ using DeviceType = typename Matrix::DeviceType;
+ using IndexType = typename Matrix::IndexType;
+ using DiagonalsShiftsType = typename Matrix::DiagonalsShiftsType;
- const IndexType rows = 8;
- const IndexType cols = 7;
+ const IndexType rows = 8;
+ const IndexType cols = 7;
- Matrix m;
- m.reset();
- m.setDimensions( rows, cols );
+ Matrix m( rows, cols );
- EXPECT_EQ( m.getRowLength( 0 ), 7 );
- EXPECT_EQ( m.getRowLength( 1 ), 7 );
- EXPECT_EQ( m.getRowLength( 2 ), 7 );
- EXPECT_EQ( m.getRowLength( 3 ), 7 );
- EXPECT_EQ( m.getRowLength( 4 ), 7 );
- EXPECT_EQ( m.getRowLength( 5 ), 7 );
- EXPECT_EQ( m.getRowLength( 6 ), 7 );
- EXPECT_EQ( m.getRowLength( 7 ), 7 );
+ EXPECT_EQ( m.getRowLength( 0 ), 2 );
+ EXPECT_EQ( m.getRowLength( 1 ), 3 );
+ EXPECT_EQ( m.getRowLength( 2 ), 3 );
+ EXPECT_EQ( m.getRowLength( 3 ), 3 );
+ EXPECT_EQ( m.getRowLength( 4 ), 3 );
+ EXPECT_EQ( m.getRowLength( 5 ), 3 );
+ EXPECT_EQ( m.getRowLength( 6 ), 2 );
+ EXPECT_EQ( m.getRowLength( 7 ), 1 );
}
template< typename Matrix >
-void test_GetNumberOfMatrixElements()
+void test_GetAllocatedElementsCount()
{
- using RealType = typename Matrix::RealType;
- using DeviceType = typename Matrix::DeviceType;
- using IndexType = typename Matrix::IndexType;
+ using RealType = typename Matrix::RealType;
+ using DeviceType = typename Matrix::DeviceType;
+ using IndexType = typename Matrix::IndexType;
- const IndexType rows = 7;
- const IndexType cols = 6;
+ const IndexType rows = 7;
+ const IndexType cols = 6;
- Matrix m;
- m.reset();
- m.setDimensions( rows, cols );
+ Matrix m( rows, cols );
- EXPECT_EQ( m.getNumberOfMatrixElements(), 42 );
+ EXPECT_EQ( m.getAllocatedElementsCount(), 21 );
}
template< typename Matrix >
void test_GetNumberOfNonzeroMatrixElements()
{
- using RealType = typename Matrix::RealType;
- using DeviceType = typename Matrix::DeviceType;
- using IndexType = typename Matrix::IndexType;
+ using RealType = typename Matrix::RealType;
+ using DeviceType = typename Matrix::DeviceType;
+ using IndexType = typename Matrix::IndexType;
-/*
- * Sets up the following 7x6 dense matrix:
- *
- * / 0 2 3 4 5 6 \
- * | 7 8 9 10 11 12 |
- * | 13 14 15 16 17 18 |
- * | 19 20 21 22 23 24 |
- * | 25 26 27 28 29 30 |
- * | 31 32 33 34 35 36 |
- * \ 37 38 39 40 41 0 /
- */
- const IndexType rows = 7;
- const IndexType cols = 6;
+ /*
+ * Sets up the following 7x6 matrix:
+ *
+ * / 0 1 0 0 0 0 \
+ * | 2 3 4 0 0 0 |
+ * | 0 5 6 7 0 0 |
+ * | 0 0 8 9 10 0 |
+ * | 0 0 0 11 12 13 |
+ * | 0 0 0 0 14 0 |
+ * \ 0 0 0 0 0 16 /
+ */
+ const IndexType rows = 7;
+ const IndexType cols = 6;
- Matrix m;
- m.reset();
- m.setDimensions( rows, cols );
+ Matrix m( rows, cols );
- RealType value = 1;
- for( IndexType i = 0; i < rows; i++ )
- for( IndexType j = 0; j < cols; j++ )
- m.setElement( i, j, value++ );
+ RealType value = 0;
+ for( IndexType i = 0; i < rows; i++ )
+ for( IndexType j = TNL::max( 0, i - 1 ); j < TNL::min( cols, i + 2 ); j++ )
+ m.setElement( i, j, value++ );
- m.setElement( 0, 0, 0); // Set the first element of the diagonal to 0.
- m.setElement( 6, 5, 0); // Set the last element of the diagonal to 0.
+ m.setElement( 5, 5, 0);
- EXPECT_EQ( m.getNumberOfNonzeroMatrixElements(), 40 );
+ EXPECT_EQ( m.getNumberOfNonzeroMatrixElements(), 15 );
}
template< typename Matrix >
void test_Reset()
{
- using RealType = typename Matrix::RealType;
- using DeviceType = typename Matrix::DeviceType;
- using IndexType = typename Matrix::IndexType;
+ using RealType = typename Matrix::RealType;
+ using DeviceType = typename Matrix::DeviceType;
+ using IndexType = typename Matrix::IndexType;
-/*
- * Sets up the following 5x4 dense matrix:
- *
- * / 0 0 0 0 \
- * | 0 0 0 0 |
- * | 0 0 0 0 |
- * | 0 0 0 0 |
- * \ 0 0 0 0 /
- */
- const IndexType rows = 5;
- const IndexType cols = 4;
+ /*
+ * Sets up the following 5x4 matrix:
+ *
+ * / 0 0 0 0 \
+ * | 0 0 0 0 |
+ * | 0 0 0 0 |
+ * | 0 0 0 0 |
+ * \ 0 0 0 0 /
+ */
+ const IndexType rows = 5;
+ const IndexType cols = 4;
- Matrix m;
- m.setDimensions( rows, cols );
+ Matrix m( rows, cols );
- m.reset();
+ m.reset();
- EXPECT_EQ( m.getRows(), 0 );
- EXPECT_EQ( m.getColumns(), 0 );
+ EXPECT_EQ( m.getRows(), 0 );
+ EXPECT_EQ( m.getColumns(), 0 );
}
template< typename Matrix >
void test_SetValue()
{
- using RealType = typename Matrix::RealType;
- using DeviceType = typename Matrix::DeviceType;
- using IndexType = typename Matrix::IndexType;
-/*
- * Sets up the following 7x6 dense matrix:
- *
- * / 1 2 3 4 5 6 \
- * | 7 8 9 10 11 12 |
- * | 13 14 15 16 17 18 |
- * | 19 20 21 22 23 24 |
- * | 25 26 27 28 29 30 |
- * | 31 32 33 34 35 36 |
- * \ 37 38 39 40 41 42 /
- */
- const IndexType rows = 7;
- const IndexType cols = 6;
+ using RealType = typename Matrix::RealType;
+ using DeviceType = typename Matrix::DeviceType;
+ using IndexType = typename Matrix::IndexType;
- Matrix m;
- m.reset();
- m.setDimensions( rows, cols );
+ /*
+ * Sets up the following 7x6 matrix:
+ *
+ * / 0 1 0 0 0 0 \
+ * | 2 3 4 0 0 0 |
+ * | 0 5 6 7 0 0 |
+ * | 0 0 8 9 10 0 |
+ * | 0 0 0 11 12 13 |
+ * | 0 0 0 0 14 0 |
+ * \ 0 0 0 0 0 16 /
+ */
+ const IndexType rows = 7;
+ const IndexType cols = 6;
- RealType value = 1;
- for( IndexType i = 0; i < rows; i++ )
- for( IndexType j = 0; j < cols; j++ )
- m.setElement( i, j, value++ );
+ Matrix m( rows, cols );
- EXPECT_EQ( m.getElement( 0, 0 ), 1 );
- EXPECT_EQ( m.getElement( 0, 1 ), 2 );
- EXPECT_EQ( m.getElement( 0, 2 ), 3 );
- EXPECT_EQ( m.getElement( 0, 3 ), 4 );
- EXPECT_EQ( m.getElement( 0, 4 ), 5 );
- EXPECT_EQ( m.getElement( 0, 5 ), 6 );
-
- EXPECT_EQ( m.getElement( 1, 0 ), 7 );
- EXPECT_EQ( m.getElement( 1, 1 ), 8 );
- EXPECT_EQ( m.getElement( 1, 2 ), 9 );
- EXPECT_EQ( m.getElement( 1, 3 ), 10 );
- EXPECT_EQ( m.getElement( 1, 4 ), 11 );
- EXPECT_EQ( m.getElement( 1, 5 ), 12 );
-
- EXPECT_EQ( m.getElement( 2, 0 ), 13 );
- EXPECT_EQ( m.getElement( 2, 1 ), 14 );
- EXPECT_EQ( m.getElement( 2, 2 ), 15 );
- EXPECT_EQ( m.getElement( 2, 3 ), 16 );
- EXPECT_EQ( m.getElement( 2, 4 ), 17 );
- EXPECT_EQ( m.getElement( 2, 5 ), 18 );
-
- EXPECT_EQ( m.getElement( 3, 0 ), 19 );
- EXPECT_EQ( m.getElement( 3, 1 ), 20 );
- EXPECT_EQ( m.getElement( 3, 2 ), 21 );
- EXPECT_EQ( m.getElement( 3, 3 ), 22 );
- EXPECT_EQ( m.getElement( 3, 4 ), 23 );
- EXPECT_EQ( m.getElement( 3, 5 ), 24 );
-
- EXPECT_EQ( m.getElement( 4, 0 ), 25 );
- EXPECT_EQ( m.getElement( 4, 1 ), 26 );
- EXPECT_EQ( m.getElement( 4, 2 ), 27 );
- EXPECT_EQ( m.getElement( 4, 3 ), 28 );
- EXPECT_EQ( m.getElement( 4, 4 ), 29 );
- EXPECT_EQ( m.getElement( 4, 5 ), 30 );
-
- EXPECT_EQ( m.getElement( 5, 0 ), 31 );
- EXPECT_EQ( m.getElement( 5, 1 ), 32 );
- EXPECT_EQ( m.getElement( 5, 2 ), 33 );
- EXPECT_EQ( m.getElement( 5, 3 ), 34 );
- EXPECT_EQ( m.getElement( 5, 4 ), 35 );
- EXPECT_EQ( m.getElement( 5, 5 ), 36 );
-
- EXPECT_EQ( m.getElement( 6, 0 ), 37 );
- EXPECT_EQ( m.getElement( 6, 1 ), 38 );
- EXPECT_EQ( m.getElement( 6, 2 ), 39 );
- EXPECT_EQ( m.getElement( 6, 3 ), 40 );
- EXPECT_EQ( m.getElement( 6, 4 ), 41 );
- EXPECT_EQ( m.getElement( 6, 5 ), 42 );
-
- // Set the values of all elements to a certain number
- m.setValue( 42 );
-
- EXPECT_EQ( m.getElement( 0, 0 ), 42 );
- EXPECT_EQ( m.getElement( 0, 1 ), 42 );
- EXPECT_EQ( m.getElement( 0, 2 ), 42 );
- EXPECT_EQ( m.getElement( 0, 3 ), 42 );
- EXPECT_EQ( m.getElement( 0, 4 ), 42 );
- EXPECT_EQ( m.getElement( 0, 5 ), 42 );
-
- EXPECT_EQ( m.getElement( 1, 0 ), 42 );
- EXPECT_EQ( m.getElement( 1, 1 ), 42 );
- EXPECT_EQ( m.getElement( 1, 2 ), 42 );
- EXPECT_EQ( m.getElement( 1, 3 ), 42 );
- EXPECT_EQ( m.getElement( 1, 4 ), 42 );
- EXPECT_EQ( m.getElement( 1, 5 ), 42 );
-
- EXPECT_EQ( m.getElement( 2, 0 ), 42 );
- EXPECT_EQ( m.getElement( 2, 1 ), 42 );
- EXPECT_EQ( m.getElement( 2, 2 ), 42 );
- EXPECT_EQ( m.getElement( 2, 3 ), 42 );
- EXPECT_EQ( m.getElement( 2, 4 ), 42 );
- EXPECT_EQ( m.getElement( 2, 5 ), 42 );
-
- EXPECT_EQ( m.getElement( 3, 0 ), 42 );
- EXPECT_EQ( m.getElement( 3, 1 ), 42 );
- EXPECT_EQ( m.getElement( 3, 2 ), 42 );
- EXPECT_EQ( m.getElement( 3, 3 ), 42 );
- EXPECT_EQ( m.getElement( 3, 4 ), 42 );
- EXPECT_EQ( m.getElement( 3, 5 ), 42 );
-
- EXPECT_EQ( m.getElement( 4, 0 ), 42 );
- EXPECT_EQ( m.getElement( 4, 1 ), 42 );
- EXPECT_EQ( m.getElement( 4, 2 ), 42 );
- EXPECT_EQ( m.getElement( 4, 3 ), 42 );
- EXPECT_EQ( m.getElement( 4, 4 ), 42 );
- EXPECT_EQ( m.getElement( 4, 5 ), 42 );
-
- EXPECT_EQ( m.getElement( 5, 0 ), 42 );
- EXPECT_EQ( m.getElement( 5, 1 ), 42 );
- EXPECT_EQ( m.getElement( 5, 2 ), 42 );
- EXPECT_EQ( m.getElement( 5, 3 ), 42 );
- EXPECT_EQ( m.getElement( 5, 4 ), 42 );
- EXPECT_EQ( m.getElement( 5, 5 ), 42 );
-
- EXPECT_EQ( m.getElement( 6, 0 ), 42 );
- EXPECT_EQ( m.getElement( 6, 1 ), 42 );
- EXPECT_EQ( m.getElement( 6, 2 ), 42 );
- EXPECT_EQ( m.getElement( 6, 3 ), 42 );
- EXPECT_EQ( m.getElement( 6, 4 ), 42 );
- EXPECT_EQ( m.getElement( 6, 5 ), 42 );
+ RealType value = 0;
+ for( IndexType i = 0; i < rows; i++ )
+ for( IndexType j = TNL::max( 0, i - 1 ); j < TNL::min( cols, i + 2 ); j++ )
+ m.setElement( i, j, value++ );
+
+ m.setElement( 5, 5, 0);
+
+ EXPECT_EQ( m.getElement( 0, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 0, 1 ), 1 );
+ EXPECT_EQ( m.getElement( 0, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 0, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 0, 4 ), 0 );
+ EXPECT_EQ( m.getElement( 0, 5 ), 0 );
+
+ EXPECT_EQ( m.getElement( 1, 0 ), 2 );
+ EXPECT_EQ( m.getElement( 1, 1 ), 3 );
+ EXPECT_EQ( m.getElement( 1, 2 ), 4 );
+ EXPECT_EQ( m.getElement( 1, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 1, 4 ), 0 );
+ EXPECT_EQ( m.getElement( 1, 5 ), 0 );
+
+ EXPECT_EQ( m.getElement( 2, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 2, 1 ), 5 );
+ EXPECT_EQ( m.getElement( 2, 2 ), 6 );
+ EXPECT_EQ( m.getElement( 2, 3 ), 7 );
+ EXPECT_EQ( m.getElement( 2, 4 ), 0 );
+ EXPECT_EQ( m.getElement( 2, 5 ), 0 );
+
+ EXPECT_EQ( m.getElement( 3, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 3, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 3, 2 ), 8 );
+ EXPECT_EQ( m.getElement( 3, 3 ), 9 );
+ EXPECT_EQ( m.getElement( 3, 4 ), 10 );
+ EXPECT_EQ( m.getElement( 3, 5 ), 0 );
+
+ EXPECT_EQ( m.getElement( 4, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 3 ), 11 );
+ EXPECT_EQ( m.getElement( 4, 4 ), 12 );
+ EXPECT_EQ( m.getElement( 4, 5 ), 13 );
+
+ EXPECT_EQ( m.getElement( 5, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 4 ), 14 );
+ EXPECT_EQ( m.getElement( 5, 5 ), 0 );
+
+ EXPECT_EQ( m.getElement( 6, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 6, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 6, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 6, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 6, 4 ), 0 );
+ EXPECT_EQ( m.getElement( 6, 5 ), 16 );
+
+ // Set the values of all elements to a certain number
+ m.setValue( 42 );
+
+ EXPECT_EQ( m.getElement( 0, 0 ), 42 );
+ EXPECT_EQ( m.getElement( 0, 1 ), 42 );
+ EXPECT_EQ( m.getElement( 0, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 0, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 0, 4 ), 0 );
+ EXPECT_EQ( m.getElement( 0, 5 ), 0 );
+
+ EXPECT_EQ( m.getElement( 1, 0 ), 42 );
+ EXPECT_EQ( m.getElement( 1, 1 ), 42 );
+ EXPECT_EQ( m.getElement( 1, 2 ), 42 );
+ EXPECT_EQ( m.getElement( 1, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 1, 4 ), 0 );
+ EXPECT_EQ( m.getElement( 1, 5 ), 0 );
+
+ EXPECT_EQ( m.getElement( 2, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 2, 1 ), 42 );
+ EXPECT_EQ( m.getElement( 2, 2 ), 42 );
+ EXPECT_EQ( m.getElement( 2, 3 ), 42 );
+ EXPECT_EQ( m.getElement( 2, 4 ), 0 );
+ EXPECT_EQ( m.getElement( 2, 5 ), 0 );
+
+ EXPECT_EQ( m.getElement( 3, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 3, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 3, 2 ), 42 );
+ EXPECT_EQ( m.getElement( 3, 3 ), 42 );
+ EXPECT_EQ( m.getElement( 3, 4 ), 42 );
+ EXPECT_EQ( m.getElement( 3, 5 ), 0 );
+
+ EXPECT_EQ( m.getElement( 4, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 3 ), 42 );
+ EXPECT_EQ( m.getElement( 4, 4 ), 42 );
+ EXPECT_EQ( m.getElement( 4, 5 ), 42 );
+
+ EXPECT_EQ( m.getElement( 5, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 4 ), 42 );
+ EXPECT_EQ( m.getElement( 5, 5 ), 42 );
+
+ EXPECT_EQ( m.getElement( 6, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 6, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 6, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 6, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 6, 4 ), 0 );
+ EXPECT_EQ( m.getElement( 6, 5 ), 42 );
}
template< typename Matrix >
void test_SetElement()
{
- using RealType = typename Matrix::RealType;
- using DeviceType = typename Matrix::DeviceType;
- using IndexType = typename Matrix::IndexType;
-/*
- * Sets up the following 5x5 dense matrix:
- *
- * / 1 2 3 4 5 \
- * | 6 7 8 9 10 |
- * | 11 12 13 14 15 |
- * | 16 17 18 19 20 |
- * \ 21 22 23 24 25 /
- */
- const IndexType rows = 5;
- const IndexType cols = 5;
+ using RealType = typename Matrix::RealType;
+ using DeviceType = typename Matrix::DeviceType;
+ using IndexType = typename Matrix::IndexType;
- Matrix m;
- m.reset();
- m.setDimensions( rows, cols );
+ /*
+ * Sets up the following 5x5 matrix:
+ *
+ * / 1 2 0 0 0 \
+ * | 6 7 8 0 0 |
+ * | 0 12 13 14 0 |
+ * | 0 0 18 19 20 |
+ * \ 0 0 0 24 25 /
+ */
+ const IndexType rows = 5;
+ const IndexType cols = 5;
- RealType value = 1;
- for( IndexType i = 0; i < rows; i++ )
- for( IndexType j = 0; j < cols; j++ )
+ Matrix m( rows, cols );
+
+ RealType value = 1;
+ for( IndexType i = 0; i < rows; i++ )
+ for( IndexType j = 0; j < cols; j++ )
+ {
+ if( abs( i - j ) > 1 )
+ {
+ EXPECT_THROW( m.setElement( i, j, value++ ), std::logic_error );
+ }
+ else
m.setElement( i, j, value++ );
+ }
+
+ EXPECT_EQ( m.getElement( 0, 0 ), 1 );
+ EXPECT_EQ( m.getElement( 0, 1 ), 2 );
+ EXPECT_EQ( m.getElement( 0, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 0, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 0, 4 ), 0 );
+
+ EXPECT_EQ( m.getElement( 1, 0 ), 6 );
+ EXPECT_EQ( m.getElement( 1, 1 ), 7 );
+ EXPECT_EQ( m.getElement( 1, 2 ), 8 );
+ EXPECT_EQ( m.getElement( 1, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 1, 4 ), 0 );
- EXPECT_EQ( m.getElement( 0, 0 ), 1 );
- EXPECT_EQ( m.getElement( 0, 1 ), 2 );
- EXPECT_EQ( m.getElement( 0, 2 ), 3 );
- EXPECT_EQ( m.getElement( 0, 3 ), 4 );
- EXPECT_EQ( m.getElement( 0, 4 ), 5 );
-
- EXPECT_EQ( m.getElement( 1, 0 ), 6 );
- EXPECT_EQ( m.getElement( 1, 1 ), 7 );
- EXPECT_EQ( m.getElement( 1, 2 ), 8 );
- EXPECT_EQ( m.getElement( 1, 3 ), 9 );
- EXPECT_EQ( m.getElement( 1, 4 ), 10 );
-
- EXPECT_EQ( m.getElement( 2, 0 ), 11 );
- EXPECT_EQ( m.getElement( 2, 1 ), 12 );
- EXPECT_EQ( m.getElement( 2, 2 ), 13 );
- EXPECT_EQ( m.getElement( 2, 3 ), 14 );
- EXPECT_EQ( m.getElement( 2, 4 ), 15 );
-
- EXPECT_EQ( m.getElement( 3, 0 ), 16 );
- EXPECT_EQ( m.getElement( 3, 1 ), 17 );
- EXPECT_EQ( m.getElement( 3, 2 ), 18 );
- EXPECT_EQ( m.getElement( 3, 3 ), 19 );
- EXPECT_EQ( m.getElement( 3, 4 ), 20 );
-
- EXPECT_EQ( m.getElement( 4, 0 ), 21 );
- EXPECT_EQ( m.getElement( 4, 1 ), 22 );
- EXPECT_EQ( m.getElement( 4, 2 ), 23 );
- EXPECT_EQ( m.getElement( 4, 3 ), 24 );
- EXPECT_EQ( m.getElement( 4, 4 ), 25 );
+ EXPECT_EQ( m.getElement( 2, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 2, 1 ), 12 );
+ EXPECT_EQ( m.getElement( 2, 2 ), 13 );
+ EXPECT_EQ( m.getElement( 2, 3 ), 14 );
+ EXPECT_EQ( m.getElement( 2, 4 ), 0 );
+
+ EXPECT_EQ( m.getElement( 3, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 3, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 3, 2 ), 18 );
+ EXPECT_EQ( m.getElement( 3, 3 ), 19 );
+ EXPECT_EQ( m.getElement( 3, 4 ), 20 );
+
+ EXPECT_EQ( m.getElement( 4, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 3 ), 24 );
+ EXPECT_EQ( m.getElement( 4, 4 ), 25 );
}
template< typename Matrix >
void test_AddElement()
{
- using RealType = typename Matrix::RealType;
- using DeviceType = typename Matrix::DeviceType;
- using IndexType = typename Matrix::IndexType;
-/*
- * Sets up the following 6x5 dense matrix:
- *
- * / 1 2 3 4 5 \
- * | 6 7 8 9 10 |
- * | 11 12 13 14 15 |
- * | 16 17 18 19 20 |
- * | 21 22 23 24 25 |
- * \ 26 27 28 29 30 /
- */
- const IndexType rows = 6;
- const IndexType cols = 5;
+ using RealType = typename Matrix::RealType;
+ using DeviceType = typename Matrix::DeviceType;
+ using IndexType = typename Matrix::IndexType;
- Matrix m;
- m.reset();
- m.setDimensions( rows, cols );
+ /*
+ * Sets up the following 6x5 matrix:
+ *
+ * / 1 2 0 0 0 \
+ * | 6 7 8 0 0 |
+ * | 0 12 13 14 0 |
+ * | 0 0 18 19 20 |
+ * | 0 0 0 24 25 |
+ * \ 0 0 0 0 30 /
+ */
+
+ const IndexType rows = 6;
+ const IndexType cols = 5;
+
+ Matrix m( rows, cols );
RealType value = 1;
for( IndexType i = 0; i < rows; i++ )
for( IndexType j = 0; j < cols; j++ )
- m.setElement( i, j, value++ );
+ {
+ if( abs( i - j ) <= 1 )
+ m.setElement( i, j, value );
+ value++;
+ }
- // Check the added elements
- EXPECT_EQ( m.getElement( 0, 0 ), 1 );
- EXPECT_EQ( m.getElement( 0, 1 ), 2 );
- EXPECT_EQ( m.getElement( 0, 2 ), 3 );
- EXPECT_EQ( m.getElement( 0, 3 ), 4 );
- EXPECT_EQ( m.getElement( 0, 4 ), 5 );
-
- EXPECT_EQ( m.getElement( 1, 0 ), 6 );
- EXPECT_EQ( m.getElement( 1, 1 ), 7 );
- EXPECT_EQ( m.getElement( 1, 2 ), 8 );
- EXPECT_EQ( m.getElement( 1, 3 ), 9 );
- EXPECT_EQ( m.getElement( 1, 4 ), 10 );
-
- EXPECT_EQ( m.getElement( 2, 0 ), 11 );
- EXPECT_EQ( m.getElement( 2, 1 ), 12 );
- EXPECT_EQ( m.getElement( 2, 2 ), 13 );
- EXPECT_EQ( m.getElement( 2, 3 ), 14 );
- EXPECT_EQ( m.getElement( 2, 4 ), 15 );
-
- EXPECT_EQ( m.getElement( 3, 0 ), 16 );
- EXPECT_EQ( m.getElement( 3, 1 ), 17 );
- EXPECT_EQ( m.getElement( 3, 2 ), 18 );
- EXPECT_EQ( m.getElement( 3, 3 ), 19 );
- EXPECT_EQ( m.getElement( 3, 4 ), 20 );
-
- EXPECT_EQ( m.getElement( 4, 0 ), 21 );
- EXPECT_EQ( m.getElement( 4, 1 ), 22 );
- EXPECT_EQ( m.getElement( 4, 2 ), 23 );
- EXPECT_EQ( m.getElement( 4, 3 ), 24 );
- EXPECT_EQ( m.getElement( 4, 4 ), 25 );
-
- EXPECT_EQ( m.getElement( 5, 0 ), 26 );
- EXPECT_EQ( m.getElement( 5, 1 ), 27 );
- EXPECT_EQ( m.getElement( 5, 2 ), 28 );
- EXPECT_EQ( m.getElement( 5, 3 ), 29 );
- EXPECT_EQ( m.getElement( 5, 4 ), 30 );
-
- // Add new elements to the old elements with a multiplying factor applied to the old elements.
-/*
- * The following setup results in the following 6x5 dense matrix:
- *
- * / 3 6 9 12 15 \
- * | 18 21 24 27 30 |
- * | 33 36 39 42 45 |
- * | 48 51 54 57 60 |
- * | 63 66 69 72 75 |
- * \ 78 81 84 87 90 /
- */
- RealType newValue = 1;
- RealType multiplicator = 2;
- for( IndexType i = 0; i < rows; i++ )
- for( IndexType j = 0; j < cols; j++ )
+ // Check the added elements
+ EXPECT_EQ( m.getElement( 0, 0 ), 1 );
+ EXPECT_EQ( m.getElement( 0, 1 ), 2 );
+ EXPECT_EQ( m.getElement( 0, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 0, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 0, 4 ), 0 );
+
+ EXPECT_EQ( m.getElement( 1, 0 ), 6 );
+ EXPECT_EQ( m.getElement( 1, 1 ), 7 );
+ EXPECT_EQ( m.getElement( 1, 2 ), 8 );
+ EXPECT_EQ( m.getElement( 1, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 1, 4 ), 0 );
+
+ EXPECT_EQ( m.getElement( 2, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 2, 1 ), 12 );
+ EXPECT_EQ( m.getElement( 2, 2 ), 13 );
+ EXPECT_EQ( m.getElement( 2, 3 ), 14 );
+ EXPECT_EQ( m.getElement( 2, 4 ), 0 );
+
+ EXPECT_EQ( m.getElement( 3, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 3, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 3, 2 ), 18 );
+ EXPECT_EQ( m.getElement( 3, 3 ), 19 );
+ EXPECT_EQ( m.getElement( 3, 4 ), 20 );
+
+ EXPECT_EQ( m.getElement( 4, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 3 ), 24 );
+ EXPECT_EQ( m.getElement( 4, 4 ), 25 );
+
+ EXPECT_EQ( m.getElement( 5, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 4 ), 30 );
+
+ // Add new elements to the old elements with a multiplying factor applied to the old elements.
+ /*
+ * The following setup results in the following 6x5 matrix:
+ *
+ * / 1 2 0 0 0 \ / 1 2 0 0 0 \ / 3 6 0 0 0 \
+ * | 6 7 8 0 0 | | 3 4 5 0 0 | | 15 18 21 0 0 |
+ * 2 * | 0 12 13 14 0 | + | 0 6 7 8 0 | = | 0 30 33 36 0 |
+ * | 0 0 18 19 20 | | 0 0 9 10 11 | | 0 0 45 48 51 |
+ * | 0 0 0 24 25 | | 0 0 0 12 13 | | 0 0 0 60 63 |
+ * \ 0 0 0 0 30 / \ 0 0 0 0 14 / \ 0 0 0 0 74 /
+ */
+
+ RealType newValue = 1;
+ RealType multiplicator = 2;
+ for( IndexType i = 0; i < rows; i++ )
+ for( IndexType j = 0; j < cols; j++ )
+ if( abs( i - j ) <= 1 )
m.addElement( i, j, newValue++, multiplicator );
- EXPECT_EQ( m.getElement( 0, 0 ), 3 );
- EXPECT_EQ( m.getElement( 0, 1 ), 6 );
- EXPECT_EQ( m.getElement( 0, 2 ), 9 );
- EXPECT_EQ( m.getElement( 0, 3 ), 12 );
- EXPECT_EQ( m.getElement( 0, 4 ), 15 );
-
- EXPECT_EQ( m.getElement( 1, 0 ), 18 );
- EXPECT_EQ( m.getElement( 1, 1 ), 21 );
- EXPECT_EQ( m.getElement( 1, 2 ), 24 );
- EXPECT_EQ( m.getElement( 1, 3 ), 27 );
- EXPECT_EQ( m.getElement( 1, 4 ), 30 );
-
- EXPECT_EQ( m.getElement( 2, 0 ), 33 );
- EXPECT_EQ( m.getElement( 2, 1 ), 36 );
- EXPECT_EQ( m.getElement( 2, 2 ), 39 );
- EXPECT_EQ( m.getElement( 2, 3 ), 42 );
- EXPECT_EQ( m.getElement( 2, 4 ), 45 );
-
- EXPECT_EQ( m.getElement( 3, 0 ), 48 );
- EXPECT_EQ( m.getElement( 3, 1 ), 51 );
- EXPECT_EQ( m.getElement( 3, 2 ), 54 );
- EXPECT_EQ( m.getElement( 3, 3 ), 57 );
- EXPECT_EQ( m.getElement( 3, 4 ), 60 );
-
- EXPECT_EQ( m.getElement( 4, 0 ), 63 );
- EXPECT_EQ( m.getElement( 4, 1 ), 66 );
- EXPECT_EQ( m.getElement( 4, 2 ), 69 );
- EXPECT_EQ( m.getElement( 4, 3 ), 72 );
- EXPECT_EQ( m.getElement( 4, 4 ), 75 );
-
- EXPECT_EQ( m.getElement( 5, 0 ), 78 );
- EXPECT_EQ( m.getElement( 5, 1 ), 81 );
- EXPECT_EQ( m.getElement( 5, 2 ), 84 );
- EXPECT_EQ( m.getElement( 5, 3 ), 87 );
- EXPECT_EQ( m.getElement( 5, 4 ), 90 );
+ EXPECT_EQ( m.getElement( 0, 0 ), 3 );
+ EXPECT_EQ( m.getElement( 0, 1 ), 6 );
+ EXPECT_EQ( m.getElement( 0, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 0, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 0, 4 ), 0 );
+
+ EXPECT_EQ( m.getElement( 1, 0 ), 15 );
+ EXPECT_EQ( m.getElement( 1, 1 ), 18 );
+ EXPECT_EQ( m.getElement( 1, 2 ), 21 );
+ EXPECT_EQ( m.getElement( 1, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 1, 4 ), 0 );
+
+ EXPECT_EQ( m.getElement( 2, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 2, 1 ), 30 );
+ EXPECT_EQ( m.getElement( 2, 2 ), 33 );
+ EXPECT_EQ( m.getElement( 2, 3 ), 36 );
+ EXPECT_EQ( m.getElement( 2, 4 ), 0 );
+
+ EXPECT_EQ( m.getElement( 3, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 3, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 3, 2 ), 45 );
+ EXPECT_EQ( m.getElement( 3, 3 ), 48 );
+ EXPECT_EQ( m.getElement( 3, 4 ), 51 );
+
+ EXPECT_EQ( m.getElement( 4, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 3 ), 60 );
+ EXPECT_EQ( m.getElement( 4, 4 ), 63 );
+
+ EXPECT_EQ( m.getElement( 5, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 4 ), 74 );
}
template< typename Matrix >
@@ -557,63 +620,56 @@ void test_SetRow()
using IndexType = typename Matrix::IndexType;
/*
- * Sets up the following 3x7 dense matrix:
+ * Sets up the following 3x7 matrix:
*
- * / 1 2 3 4 5 6 7 \
- * | 8 9 10 11 12 13 14 |
- * \ 15 16 17 18 19 20 21 /
+ * / 1 2 0 0 0 0 0 \
+ * | 8 9 10 0 0 0 0 |
+ * \ 0 16 17 18 0 0 0 /
*/
const IndexType rows = 3;
const IndexType cols = 7;
- Matrix m;
- m.reset();
- m.setDimensions( rows, cols );
-
- RealType value = 1;
- for( IndexType i = 0; i < rows; i++ )
- for( IndexType j = 0; j < cols; j++ )
- m.setElement( i, j, value++ );
+ Matrix m( rows, cols );
auto matrix_view = m.getView();
auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
- RealType values[ 3 ][ 5 ] {
- { 11, 11, 11, 11, 11 },
- { 22, 22, 22, 22, 22 },
- { 33, 33, 33, 33, 33 } };
- IndexType columnIndexes[ 3 ][ 5 ] {
- { 0, 1, 2, 3, 4 },
- { 0, 1, 2, 3, 4 },
- { 2, 3, 4, 5, 6 } };
+ RealType values[ 3 ][ 3 ] {
+ { 1, 2, 0 },
+ { 8, 9, 10 },
+ { 16, 17, 18 } };
auto row = matrix_view.getRow( rowIdx );
- for( IndexType i = 0; i < 5; i++ )
- row.setElement( i, values[ rowIdx ][ i ] );
+ for( IndexType i = 0; i < 3; i++ )
+ {
+ if( rowIdx == 0 && i > 1 )
+ break;
+ row.setElement( i, values[ rowIdx ][ i ] );
+ }
};
TNL::Algorithms::ParallelFor< DeviceType >::exec( 0, 3, f );
- EXPECT_EQ( m.getElement( 0, 0 ), 11 );
- EXPECT_EQ( m.getElement( 0, 1 ), 11 );
- EXPECT_EQ( m.getElement( 0, 2 ), 11 );
- EXPECT_EQ( m.getElement( 0, 3 ), 11 );
- EXPECT_EQ( m.getElement( 0, 4 ), 11 );
- EXPECT_EQ( m.getElement( 0, 5 ), 6 );
- EXPECT_EQ( m.getElement( 0, 6 ), 7 );
-
- EXPECT_EQ( m.getElement( 1, 0 ), 22 );
- EXPECT_EQ( m.getElement( 1, 1 ), 22 );
- EXPECT_EQ( m.getElement( 1, 2 ), 22 );
- EXPECT_EQ( m.getElement( 1, 3 ), 22 );
- EXPECT_EQ( m.getElement( 1, 4 ), 22 );
- EXPECT_EQ( m.getElement( 1, 5 ), 13 );
- EXPECT_EQ( m.getElement( 1, 6 ), 14 );
-
- EXPECT_EQ( m.getElement( 2, 0 ), 15 );
+ EXPECT_EQ( m.getElement( 0, 0 ), 1 );
+ EXPECT_EQ( m.getElement( 0, 1 ), 2 );
+ EXPECT_EQ( m.getElement( 0, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 0, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 0, 4 ), 0 );
+ EXPECT_EQ( m.getElement( 0, 5 ), 0 );
+ EXPECT_EQ( m.getElement( 0, 6 ), 0 );
+
+ EXPECT_EQ( m.getElement( 1, 0 ), 8 );
+ EXPECT_EQ( m.getElement( 1, 1 ), 9 );
+ EXPECT_EQ( m.getElement( 1, 2 ), 10 );
+ EXPECT_EQ( m.getElement( 1, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 1, 4 ), 0 );
+ EXPECT_EQ( m.getElement( 1, 5 ), 0 );
+ EXPECT_EQ( m.getElement( 1, 6 ), 0 );
+
+ EXPECT_EQ( m.getElement( 2, 0 ), 0 );
EXPECT_EQ( m.getElement( 2, 1 ), 16 );
- EXPECT_EQ( m.getElement( 2, 2 ), 33 );
- EXPECT_EQ( m.getElement( 2, 3 ), 33 );
- EXPECT_EQ( m.getElement( 2, 4 ), 33 );
- EXPECT_EQ( m.getElement( 2, 5 ), 33 );
- EXPECT_EQ( m.getElement( 2, 6 ), 33 );
+ EXPECT_EQ( m.getElement( 2, 2 ), 17 );
+ EXPECT_EQ( m.getElement( 2, 3 ), 18 );
+ EXPECT_EQ( m.getElement( 2, 4 ), 0 );
+ EXPECT_EQ( m.getElement( 2, 5 ), 0 );
+ EXPECT_EQ( m.getElement( 2, 6 ), 0 );
}
template< typename Matrix >
@@ -623,14 +679,14 @@ void test_AddRow()
using DeviceType = typename Matrix::DeviceType;
using IndexType = typename Matrix::IndexType;
/*
- * Sets up the following 6x5 dense matrix:
+ * Sets up the following 6x5 matrix:
*
- * / 1 2 3 4 5 \
- * | 6 7 8 9 10 |
- * | 11 12 13 14 15 |
- * | 16 17 18 19 20 |
- * | 21 22 23 24 25 |
- * \ 26 27 28 29 30 /
+ * / 1 2 0 0 0 \
+ * | 6 7 8 0 0 |
+ * | 0 12 13 14 0 |
+ * | 0 0 18 19 20 |
+ * | 0 0 0 24 25 |
+ * \ 0 0 0 0 30 /
*/
const IndexType rows = 6;
@@ -641,68 +697,72 @@ void test_AddRow()
RealType value = 1;
for( IndexType i = 0; i < rows; i++ )
for( IndexType j = 0; j < cols; j++ )
- m.setElement( i, j, value++ );
+ {
+ if( abs( i - j ) <= 1 )
+ m.setElement( i, j, value );
+ value++;
+ }
// Check the added elements
EXPECT_EQ( m.getElement( 0, 0 ), 1 );
EXPECT_EQ( m.getElement( 0, 1 ), 2 );
- EXPECT_EQ( m.getElement( 0, 2 ), 3 );
- EXPECT_EQ( m.getElement( 0, 3 ), 4 );
- EXPECT_EQ( m.getElement( 0, 4 ), 5 );
+ EXPECT_EQ( m.getElement( 0, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 0, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 0, 4 ), 0 );
EXPECT_EQ( m.getElement( 1, 0 ), 6 );
EXPECT_EQ( m.getElement( 1, 1 ), 7 );
EXPECT_EQ( m.getElement( 1, 2 ), 8 );
- EXPECT_EQ( m.getElement( 1, 3 ), 9 );
- EXPECT_EQ( m.getElement( 1, 4 ), 10 );
+ EXPECT_EQ( m.getElement( 1, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 1, 4 ), 0 );
- EXPECT_EQ( m.getElement( 2, 0 ), 11 );
+ EXPECT_EQ( m.getElement( 2, 0 ), 0 );
EXPECT_EQ( m.getElement( 2, 1 ), 12 );
EXPECT_EQ( m.getElement( 2, 2 ), 13 );
EXPECT_EQ( m.getElement( 2, 3 ), 14 );
- EXPECT_EQ( m.getElement( 2, 4 ), 15 );
+ EXPECT_EQ( m.getElement( 2, 4 ), 0 );
- EXPECT_EQ( m.getElement( 3, 0 ), 16 );
- EXPECT_EQ( m.getElement( 3, 1 ), 17 );
+ EXPECT_EQ( m.getElement( 3, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 3, 1 ), 0 );
EXPECT_EQ( m.getElement( 3, 2 ), 18 );
EXPECT_EQ( m.getElement( 3, 3 ), 19 );
EXPECT_EQ( m.getElement( 3, 4 ), 20 );
- EXPECT_EQ( m.getElement( 4, 0 ), 21 );
- EXPECT_EQ( m.getElement( 4, 1 ), 22 );
- EXPECT_EQ( m.getElement( 4, 2 ), 23 );
+ EXPECT_EQ( m.getElement( 4, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 2 ), 0 );
EXPECT_EQ( m.getElement( 4, 3 ), 24 );
EXPECT_EQ( m.getElement( 4, 4 ), 25 );
- EXPECT_EQ( m.getElement( 5, 0 ), 26 );
- EXPECT_EQ( m.getElement( 5, 1 ), 27 );
- EXPECT_EQ( m.getElement( 5, 2 ), 28 );
- EXPECT_EQ( m.getElement( 5, 3 ), 29 );
+ EXPECT_EQ( m.getElement( 5, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 3 ), 0 );
EXPECT_EQ( m.getElement( 5, 4 ), 30 );
// Add new elements to the old elements with a multiplying factor applied to the old elements.
/*
* The following setup results in the following 6x5 sparse matrix:
*
- * / 3 6 9 12 15 \
- * | 18 21 24 27 30 |
- * | 33 36 39 42 45 |
- * | 48 51 54 57 60 |
- * | 63 66 69 72 75 |
- * \ 78 81 84 87 90 /
+ * / 0 0 0 0 0 0 \ / 1 2 0 0 0 \ / 11 11 0 0 0 \ / 11 11 0 0 0 \
+ * | 0 1 0 0 0 0 | | 6 7 8 0 0 | | 22 22 22 0 0 | | 28 29 30 0 0 |
+ * | 0 0 2 0 0 0 | * | 0 12 13 14 0 | + | 0 33 33 33 0 | = | 0 57 59 61 0 |
+ * | 0 0 0 3 0 0 | | 0 0 18 19 20 | | 0 0 44 44 44 | | 0 0 98 101 104 |
+ * | 0 0 0 0 4 0 | | 0 0 0 24 25 | | 0 0 0 55 55 | | 0 0 0 151 155 |
+ * \ 0 0 0 0 0 5 / \ 0 0 0 0 30 / \ 0 0 0 0 66 / \ 0 0 0 0 216 /
*/
auto matrix_view = m.getView();
auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
- RealType values[ 6 ][ 5 ] {
- { 11, 11, 11, 11, 0 },
- { 22, 22, 22, 22, 0 },
- { 33, 33, 33, 33, 0 },
- { 44, 44, 44, 44, 0 },
- { 55, 55, 55, 55, 0 },
- { 66, 66, 66, 66, 0 } };
+ RealType values[ 6 ][ 3 ] {
+ { 11, 11, 0 },
+ { 22, 22, 22 },
+ { 33, 33, 33 },
+ { 44, 44, 44 },
+ { 55, 55, 55 },
+ { 66, 66, 66 } };
auto row = matrix_view.getRow( rowIdx );
- for( IndexType i = 0; i < 5; i++ )
+ for( IndexType i = 0; i < 3; i++ )
{
RealType& val = row.getValue( i );
val = rowIdx * val + values[ rowIdx ][ i ];
@@ -711,208 +771,207 @@ void test_AddRow()
TNL::Algorithms::ParallelFor< DeviceType >::exec( 0, 6, f );
- EXPECT_EQ( m.getElement( 0, 0 ), 11 );
- EXPECT_EQ( m.getElement( 0, 1 ), 11 );
- EXPECT_EQ( m.getElement( 0, 2 ), 11 );
- EXPECT_EQ( m.getElement( 0, 3 ), 11 );
- EXPECT_EQ( m.getElement( 0, 4 ), 0 );
-
- EXPECT_EQ( m.getElement( 1, 0 ), 28 );
- EXPECT_EQ( m.getElement( 1, 1 ), 29 );
- EXPECT_EQ( m.getElement( 1, 2 ), 30 );
- EXPECT_EQ( m.getElement( 1, 3 ), 31 );
- EXPECT_EQ( m.getElement( 1, 4 ), 10 );
-
- EXPECT_EQ( m.getElement( 2, 0 ), 55 );
- EXPECT_EQ( m.getElement( 2, 1 ), 57 );
- EXPECT_EQ( m.getElement( 2, 2 ), 59 );
- EXPECT_EQ( m.getElement( 2, 3 ), 61 );
- EXPECT_EQ( m.getElement( 2, 4 ), 30 );
-
- EXPECT_EQ( m.getElement( 3, 0 ), 92 );
- EXPECT_EQ( m.getElement( 3, 1 ), 95 );
- EXPECT_EQ( m.getElement( 3, 2 ), 98 );
- EXPECT_EQ( m.getElement( 3, 3 ), 101 );
- EXPECT_EQ( m.getElement( 3, 4 ), 60 );
-
- EXPECT_EQ( m.getElement( 4, 0 ), 139 );
- EXPECT_EQ( m.getElement( 4, 1 ), 143 );
- EXPECT_EQ( m.getElement( 4, 2 ), 147 );
- EXPECT_EQ( m.getElement( 4, 3 ), 151 );
- EXPECT_EQ( m.getElement( 4, 4 ), 100 );
-
- EXPECT_EQ( m.getElement( 5, 0 ), 196 );
- EXPECT_EQ( m.getElement( 5, 1 ), 201 );
- EXPECT_EQ( m.getElement( 5, 2 ), 206 );
- EXPECT_EQ( m.getElement( 5, 3 ), 211 );
- EXPECT_EQ( m.getElement( 5, 4 ), 150 );
+ EXPECT_EQ( m.getElement( 0, 0 ), 11 );
+ EXPECT_EQ( m.getElement( 0, 1 ), 11 );
+ EXPECT_EQ( m.getElement( 0, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 0, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 0, 4 ), 0 );
+
+ EXPECT_EQ( m.getElement( 1, 0 ), 28 );
+ EXPECT_EQ( m.getElement( 1, 1 ), 29 );
+ EXPECT_EQ( m.getElement( 1, 2 ), 30 );
+ EXPECT_EQ( m.getElement( 1, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 1, 4 ), 0 );
+
+ EXPECT_EQ( m.getElement( 2, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 2, 1 ), 57 );
+ EXPECT_EQ( m.getElement( 2, 2 ), 59 );
+ EXPECT_EQ( m.getElement( 2, 3 ), 61 );
+ EXPECT_EQ( m.getElement( 2, 4 ), 0 );
+
+ EXPECT_EQ( m.getElement( 3, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 3, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 3, 2 ), 98 );
+ EXPECT_EQ( m.getElement( 3, 3 ), 101 );
+ EXPECT_EQ( m.getElement( 3, 4 ), 104 );
+
+ EXPECT_EQ( m.getElement( 4, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 3 ), 151 );
+ EXPECT_EQ( m.getElement( 4, 4 ), 155 );
+
+ EXPECT_EQ( m.getElement( 5, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 4 ), 216 );
}
template< typename Matrix >
void test_VectorProduct()
{
- using RealType = typename Matrix::RealType;
- using DeviceType = typename Matrix::DeviceType;
- using IndexType = typename Matrix::IndexType;
-/*
- * Sets up the following 5x4 dense matrix:
- *
- * / 1 2 3 4 \
- * | 5 6 7 8 |
- * | 9 10 11 12 |
- * | 13 14 15 16 |
- * \ 17 18 19 20 /
- */
- const IndexType rows = 5;
- const IndexType cols = 4;
+ using RealType = typename Matrix::RealType;
+ using DeviceType = typename Matrix::DeviceType;
+ using IndexType = typename Matrix::IndexType;
- Matrix m;
- m.reset();
- m.setDimensions( rows, cols );
+ /*
+ * Sets up the following 5x4 matrix:
+ *
+ * / 1 2 0 0 \
+ * | 5 6 7 0 |
+ * | 0 10 11 12 |
+ * | 0 0 15 16 |
+ * \ 0 0 0 20 /
+ */
+ const IndexType rows = 5;
+ const IndexType cols = 4;
- RealType value = 1;
- for( IndexType i = 0; i < rows; i++ )
- for( IndexType j = 0; j < cols; j++)
- m.setElement( i, j, value++ );
+ Matrix m( rows, cols );
- using VectorType = TNL::Containers::Vector< RealType, DeviceType, IndexType >;
+ RealType value = 1;
+ for( IndexType i = 0; i < rows; i++ )
+ for( IndexType j = 0; j < cols; j++)
+ {
+ if( abs( i - j ) <= 1 )
+ m.setElement( i, j, value );
+ value++;
+ }
- VectorType inVector;
- inVector.setSize( 4 );
- for( IndexType i = 0; i < inVector.getSize(); i++ )
- inVector.setElement( i, 2 );
+ using VectorType = TNL::Containers::Vector< RealType, DeviceType, IndexType >;
- VectorType outVector;
- outVector.setSize( 5 );
- for( IndexType j = 0; j < outVector.getSize(); j++ )
- outVector.setElement( j, 0 );
+ VectorType inVector( 4 );
+ inVector = 2;
+ VectorType outVector( 5 );
+ outVector = 0;
- m.vectorProduct( inVector, outVector);
+ m.vectorProduct( inVector, outVector);
- EXPECT_EQ( outVector.getElement( 0 ), 20 );
- EXPECT_EQ( outVector.getElement( 1 ), 52 );
- EXPECT_EQ( outVector.getElement( 2 ), 84 );
- EXPECT_EQ( outVector.getElement( 3 ), 116 );
- EXPECT_EQ( outVector.getElement( 4 ), 148 );
+ EXPECT_EQ( outVector.getElement( 0 ), 6 );
+ EXPECT_EQ( outVector.getElement( 1 ), 36 );
+ EXPECT_EQ( outVector.getElement( 2 ), 66 );
+ EXPECT_EQ( outVector.getElement( 3 ), 62 );
+ EXPECT_EQ( outVector.getElement( 4 ), 40 );
}
-template< typename Matrix >
+template< typename Matrix1, typename Matrix2 = Matrix1 >
void test_AddMatrix()
{
- using RealType = typename Matrix::RealType;
- using DeviceType = typename Matrix::DeviceType;
- using IndexType = typename Matrix::IndexType;
-/*
- * Sets up the following 5x4 dense matrix:
- *
- * / 1 2 3 4 \
- * | 5 6 7 8 |
- * | 9 10 11 12 |
- * | 13 14 15 16 |
- * \ 17 18 19 20 /
- */
- const IndexType rows = 5;
- const IndexType cols = 4;
+ using RealType = typename Matrix1::RealType;
+ using DeviceType = typename Matrix1::DeviceType;
+ using IndexType = typename Matrix1::IndexType;
- Matrix m;
- m.reset();
- m.setDimensions( rows, cols );
+ /*
+ * Sets up the following 5x4 matrix:
+ *
+ * / 1 2 0 0 \
+ * | 5 6 7 0 |
+ * | 0 10 11 12 |
+ * | 0 0 15 16 |
+ * \ 0 0 0 20 /
+ */
+ const IndexType rows = 5;
+ const IndexType cols = 4;
- RealType value = 1;
- for( IndexType i = 0; i < rows; i++ )
- for( IndexType j = 0; j < cols; j++)
- m.setElement( i, j, value++ );
+ Matrix1 m( rows, cols );
-/*
- * Sets up the following 5x4 dense matrix:
- *
- * / 1 2 3 4 \
- * | 5 6 7 8 |
- * | 9 10 11 12 |
- * | 13 14 15 16 |
- * \ 17 18 19 20 /
- */
+ RealType value = 1;
+ for( IndexType i = 0; i < rows; i++ )
+ for( IndexType j = 0; j < cols; j++)
+ {
+ if( abs( i - j ) <= 1 )
+ m.setElement( i, j, value );
+ value++;
+ }
- Matrix m2;
- m2.reset();
- m2.setDimensions( rows, cols );
+ /*
+ * Sets up the following 5x4 matrix:
+ *
+ * / 1 2 0 0 \
+ * | 3 4 5 0 |
+ * | 0 6 7 8 |
+ * | 0 0 9 10 |
+ * \ 0 0 0 11 /
+ */
+ Matrix2 m2( rows, cols );
- RealType newValue = 1;
- for( IndexType i = 0; i < rows; i++ )
- for( IndexType j = 0; j < cols; j++)
+ RealType newValue = 1;
+ for( IndexType i = 0; i < rows; i++ )
+ for( IndexType j = 0; j < cols; j++)
+ if( abs( i - j ) <= 1 )
m2.setElement( i, j, newValue++ );
- /*
- * Sets up the following 5x4 dense matrix:
- *
- * / 1 2 3 4 \
- * | 5 6 7 8 |
- * | 9 10 11 12 |
- * | 13 14 15 16 |
- * \ 17 18 19 20 /
- */
+ /*
+ * Compute the following 5x4 matrix:
+ *
+ * / 1 2 0 0 \ / 1 2 0 0 \ / 3 6 0 0 \
+ * | 5 6 7 0 | | 3 4 5 0 | | 11 14 17 0 |
+ * | 0 10 11 12 | + 2 * | 0 6 7 8 | = | 0 22 25 28 |
+ * | 0 0 15 16 | | 0 0 9 10 | | 0 0 33 36 |
+ * \ 0 0 0 20 / \ 0 0 0 11 / \ 0 0 0 42 /
+ */
- Matrix mResult;
- mResult.reset();
- mResult.setDimensions( rows, cols );
-
- mResult = m;
-
- RealType matrixMultiplicator = 2;
- RealType thisMatrixMultiplicator = 1;
-
- mResult.addMatrix( m2, matrixMultiplicator, thisMatrixMultiplicator );
-
- EXPECT_EQ( mResult.getElement( 0, 0 ), matrixMultiplicator * m2.getElement( 0, 0 ) + thisMatrixMultiplicator * m.getElement( 0, 0 ) );
- EXPECT_EQ( mResult.getElement( 0, 1 ), matrixMultiplicator * m2.getElement( 0, 1 ) + thisMatrixMultiplicator * m.getElement( 0, 1 ) );
- EXPECT_EQ( mResult.getElement( 0, 2 ), matrixMultiplicator * m2.getElement( 0, 2 ) + thisMatrixMultiplicator * m.getElement( 0, 2 ) );
- EXPECT_EQ( mResult.getElement( 0, 3 ), matrixMultiplicator * m2.getElement( 0, 3 ) + thisMatrixMultiplicator * m.getElement( 0, 3 ) );
-
- EXPECT_EQ( mResult.getElement( 1, 0 ), matrixMultiplicator * m2.getElement( 1, 0 ) + thisMatrixMultiplicator * m.getElement( 1, 0 ) );
- EXPECT_EQ( mResult.getElement( 1, 1 ), matrixMultiplicator * m2.getElement( 1, 1 ) + thisMatrixMultiplicator * m.getElement( 1, 1 ) );
- EXPECT_EQ( mResult.getElement( 1, 2 ), matrixMultiplicator * m2.getElement( 1, 2 ) + thisMatrixMultiplicator * m.getElement( 1, 2 ) );
- EXPECT_EQ( mResult.getElement( 1, 3 ), matrixMultiplicator * m2.getElement( 1, 3 ) + thisMatrixMultiplicator * m.getElement( 1, 3 ) );
-
- EXPECT_EQ( mResult.getElement( 2, 0 ), matrixMultiplicator * m2.getElement( 2, 0 ) + thisMatrixMultiplicator * m.getElement( 2, 0 ) );
- EXPECT_EQ( mResult.getElement( 2, 1 ), matrixMultiplicator * m2.getElement( 2, 1 ) + thisMatrixMultiplicator * m.getElement( 2, 1 ) );
- EXPECT_EQ( mResult.getElement( 2, 2 ), matrixMultiplicator * m2.getElement( 2, 2 ) + thisMatrixMultiplicator * m.getElement( 2, 2 ) );
- EXPECT_EQ( mResult.getElement( 2, 3 ), matrixMultiplicator * m2.getElement( 2, 3 ) + thisMatrixMultiplicator * m.getElement( 2, 3 ) );
-
- EXPECT_EQ( mResult.getElement( 3, 0 ), matrixMultiplicator * m2.getElement( 3, 0 ) + thisMatrixMultiplicator * m.getElement( 3, 0 ) );
- EXPECT_EQ( mResult.getElement( 3, 1 ), matrixMultiplicator * m2.getElement( 3, 1 ) + thisMatrixMultiplicator * m.getElement( 3, 1 ) );
- EXPECT_EQ( mResult.getElement( 3, 2 ), matrixMultiplicator * m2.getElement( 3, 2 ) + thisMatrixMultiplicator * m.getElement( 3, 2 ) );
- EXPECT_EQ( mResult.getElement( 3, 3 ), matrixMultiplicator * m2.getElement( 3, 3 ) + thisMatrixMultiplicator * m.getElement( 3, 3 ) );
-
- EXPECT_EQ( mResult.getElement( 4, 0 ), matrixMultiplicator * m2.getElement( 4, 0 ) + thisMatrixMultiplicator * m.getElement( 4, 0 ) );
- EXPECT_EQ( mResult.getElement( 4, 1 ), matrixMultiplicator * m2.getElement( 4, 1 ) + thisMatrixMultiplicator * m.getElement( 4, 1 ) );
- EXPECT_EQ( mResult.getElement( 4, 2 ), matrixMultiplicator * m2.getElement( 4, 2 ) + thisMatrixMultiplicator * m.getElement( 4, 2 ) );
- EXPECT_EQ( mResult.getElement( 4, 3 ), matrixMultiplicator * m2.getElement( 4, 3 ) + thisMatrixMultiplicator * m.getElement( 4, 3 ) );
-
- EXPECT_EQ( mResult.getElement( 0, 0 ), 3 );
- EXPECT_EQ( mResult.getElement( 0, 1 ), 6 );
- EXPECT_EQ( mResult.getElement( 0, 2 ), 9 );
- EXPECT_EQ( mResult.getElement( 0, 3 ), 12 );
-
- EXPECT_EQ( mResult.getElement( 1, 0 ), 15 );
- EXPECT_EQ( mResult.getElement( 1, 1 ), 18 );
- EXPECT_EQ( mResult.getElement( 1, 2 ), 21 );
- EXPECT_EQ( mResult.getElement( 1, 3 ), 24 );
-
- EXPECT_EQ( mResult.getElement( 2, 0 ), 27 );
- EXPECT_EQ( mResult.getElement( 2, 1 ), 30 );
- EXPECT_EQ( mResult.getElement( 2, 2 ), 33 );
- EXPECT_EQ( mResult.getElement( 2, 3 ), 36 );
-
- EXPECT_EQ( mResult.getElement( 3, 0 ), 39 );
- EXPECT_EQ( mResult.getElement( 3, 1 ), 42 );
- EXPECT_EQ( mResult.getElement( 3, 2 ), 45 );
- EXPECT_EQ( mResult.getElement( 3, 3 ), 48 );
-
- EXPECT_EQ( mResult.getElement( 4, 0 ), 51 );
- EXPECT_EQ( mResult.getElement( 4, 1 ), 54 );
- EXPECT_EQ( mResult.getElement( 4, 2 ), 57 );
- EXPECT_EQ( mResult.getElement( 4, 3 ), 60 );
+ Matrix1 mResult;
+ mResult.reset();
+ mResult.setDimensions( rows, cols );
+
+ mResult = m;
+
+ RealType matrixMultiplicator = 2;
+ RealType thisMatrixMultiplicator = 1;
+
+ mResult.addMatrix( m2, matrixMultiplicator, thisMatrixMultiplicator );
+
+ EXPECT_EQ( mResult.getElement( 0, 0 ), matrixMultiplicator * m2.getElement( 0, 0 ) + thisMatrixMultiplicator * m.getElement( 0, 0 ) );
+ EXPECT_EQ( mResult.getElement( 0, 1 ), matrixMultiplicator * m2.getElement( 0, 1 ) + thisMatrixMultiplicator * m.getElement( 0, 1 ) );
+ EXPECT_EQ( mResult.getElement( 0, 2 ), matrixMultiplicator * m2.getElement( 0, 2 ) + thisMatrixMultiplicator * m.getElement( 0, 2 ) );
+ EXPECT_EQ( mResult.getElement( 0, 3 ), matrixMultiplicator * m2.getElement( 0, 3 ) + thisMatrixMultiplicator * m.getElement( 0, 3 ) );
+
+ EXPECT_EQ( mResult.getElement( 1, 0 ), matrixMultiplicator * m2.getElement( 1, 0 ) + thisMatrixMultiplicator * m.getElement( 1, 0 ) );
+ EXPECT_EQ( mResult.getElement( 1, 1 ), matrixMultiplicator * m2.getElement( 1, 1 ) + thisMatrixMultiplicator * m.getElement( 1, 1 ) );
+ EXPECT_EQ( mResult.getElement( 1, 2 ), matrixMultiplicator * m2.getElement( 1, 2 ) + thisMatrixMultiplicator * m.getElement( 1, 2 ) );
+ EXPECT_EQ( mResult.getElement( 1, 3 ), matrixMultiplicator * m2.getElement( 1, 3 ) + thisMatrixMultiplicator * m.getElement( 1, 3 ) );
+
+ EXPECT_EQ( mResult.getElement( 2, 0 ), matrixMultiplicator * m2.getElement( 2, 0 ) + thisMatrixMultiplicator * m.getElement( 2, 0 ) );
+ EXPECT_EQ( mResult.getElement( 2, 1 ), matrixMultiplicator * m2.getElement( 2, 1 ) + thisMatrixMultiplicator * m.getElement( 2, 1 ) );
+ EXPECT_EQ( mResult.getElement( 2, 2 ), matrixMultiplicator * m2.getElement( 2, 2 ) + thisMatrixMultiplicator * m.getElement( 2, 2 ) );
+ EXPECT_EQ( mResult.getElement( 2, 3 ), matrixMultiplicator * m2.getElement( 2, 3 ) + thisMatrixMultiplicator * m.getElement( 2, 3 ) );
+
+ EXPECT_EQ( mResult.getElement( 3, 0 ), matrixMultiplicator * m2.getElement( 3, 0 ) + thisMatrixMultiplicator * m.getElement( 3, 0 ) );
+ EXPECT_EQ( mResult.getElement( 3, 1 ), matrixMultiplicator * m2.getElement( 3, 1 ) + thisMatrixMultiplicator * m.getElement( 3, 1 ) );
+ EXPECT_EQ( mResult.getElement( 3, 2 ), matrixMultiplicator * m2.getElement( 3, 2 ) + thisMatrixMultiplicator * m.getElement( 3, 2 ) );
+ EXPECT_EQ( mResult.getElement( 3, 3 ), matrixMultiplicator * m2.getElement( 3, 3 ) + thisMatrixMultiplicator * m.getElement( 3, 3 ) );
+
+ EXPECT_EQ( mResult.getElement( 4, 0 ), matrixMultiplicator * m2.getElement( 4, 0 ) + thisMatrixMultiplicator * m.getElement( 4, 0 ) );
+ EXPECT_EQ( mResult.getElement( 4, 1 ), matrixMultiplicator * m2.getElement( 4, 1 ) + thisMatrixMultiplicator * m.getElement( 4, 1 ) );
+ EXPECT_EQ( mResult.getElement( 4, 2 ), matrixMultiplicator * m2.getElement( 4, 2 ) + thisMatrixMultiplicator * m.getElement( 4, 2 ) );
+ EXPECT_EQ( mResult.getElement( 4, 3 ), matrixMultiplicator * m2.getElement( 4, 3 ) + thisMatrixMultiplicator * m.getElement( 4, 3 ) );
+
+ EXPECT_EQ( mResult.getElement( 0, 0 ), 3 );
+ EXPECT_EQ( mResult.getElement( 0, 1 ), 6 );
+ EXPECT_EQ( mResult.getElement( 0, 2 ), 0 );
+ EXPECT_EQ( mResult.getElement( 0, 3 ), 0 );
+
+ EXPECT_EQ( mResult.getElement( 1, 0 ), 11 );
+ EXPECT_EQ( mResult.getElement( 1, 1 ), 14 );
+ EXPECT_EQ( mResult.getElement( 1, 2 ), 17 );
+ EXPECT_EQ( mResult.getElement( 1, 3 ), 0 );
+
+ EXPECT_EQ( mResult.getElement( 2, 0 ), 0 );
+ EXPECT_EQ( mResult.getElement( 2, 1 ), 22 );
+ EXPECT_EQ( mResult.getElement( 2, 2 ), 25 );
+ EXPECT_EQ( mResult.getElement( 2, 3 ), 28 );
+
+ EXPECT_EQ( mResult.getElement( 3, 0 ), 0 );
+ EXPECT_EQ( mResult.getElement( 3, 1 ), 0 );
+ EXPECT_EQ( mResult.getElement( 3, 2 ), 33 );
+ EXPECT_EQ( mResult.getElement( 3, 3 ), 36 );
+
+ EXPECT_EQ( mResult.getElement( 4, 0 ), 0 );
+ EXPECT_EQ( mResult.getElement( 4, 1 ), 0 );
+ EXPECT_EQ( mResult.getElement( 4, 2 ), 0 );
+ EXPECT_EQ( mResult.getElement( 4, 3 ), 42 );
}
template< typename Matrix >
@@ -922,7 +981,7 @@ void test_GetMatrixProduct()
using DeviceType = typename Matrix::DeviceType;
using IndexType = typename Matrix::IndexType;
/*
- * Sets up the following 5x4 dense matrix:
+ * Sets up the following 5x4 matrix:
*
* / 1 2 3 4 \
* | 5 6 7 8 |
@@ -943,7 +1002,7 @@ void test_GetMatrixProduct()
leftMatrix.setElement( i, j, value++ );
/*
- * Sets up the following 4x5 dense matrix:
+ * Sets up the following 4x5 matrix:
*
* / 1 2 3 4 5 \
* | 6 7 8 9 10 |
@@ -963,7 +1022,7 @@ void test_GetMatrixProduct()
rightMatrix.setElement( i, j, newValue++ );
/*
- * Sets up the following 5x5 resulting dense matrix:
+ * Sets up the following 5x5 resulting matrix:
*
* / 0 0 0 0 \
* | 0 0 0 0 |
@@ -1027,7 +1086,7 @@ void test_GetTransposition()
using DeviceType = typename Matrix::DeviceType;
using IndexType = typename Matrix::IndexType;
/*
- * Sets up the following 3x2 dense matrix:
+ * Sets up the following 3x2 matrix:
*
* / 1 2 \
* | 3 4 |
@@ -1048,7 +1107,7 @@ void test_GetTransposition()
m.print( std::cout );
/*
- * Sets up the following 2x3 dense matrix:
+ * Sets up the following 2x3 matrix:
*
* / 0 0 0 \
* \ 0 0 0 /
@@ -1066,7 +1125,7 @@ void test_GetTransposition()
mTransposed.print( std::cout );
/*
- * Should result in the following 2x3 dense matrix:
+ * Should result in the following 2x3 matrix:
*
* / 1 3 5 \
* \ 2 4 6 /
@@ -1089,7 +1148,7 @@ void test_PerformSORIteration()
using DeviceType = typename Matrix::DeviceType;
using IndexType = typename Matrix::IndexType;
/*
- * Sets up the following 4x4 dense matrix:
+ * Sets up the following 4x4 matrix:
*
* / 4 1 1 1 \
* | 1 4 1 1 |
@@ -1164,43 +1223,44 @@ void test_AssignmentOperator()
using RealType = typename Matrix::RealType;
using DeviceType = typename Matrix::DeviceType;
using IndexType = typename Matrix::IndexType;
+ constexpr bool rowMajorOrder = Matrix::getRowMajorOrder();
- using MultidiagonalHost = TNL::Matrices::Multidiagonal< RealType, TNL::Devices::Host, IndexType >;
- using MultidiagonalCuda = TNL::Matrices::Multidiagonal< RealType, TNL::Devices::Cuda, IndexType >;
+ using MultidiagonalHost = TNL::Matrices::Multidiagonal< RealType, TNL::Devices::Host, IndexType, rowMajorOrder >;
+ using MultidiagonalCuda = TNL::Matrices::Multidiagonal< RealType, TNL::Devices::Cuda, IndexType, !rowMajorOrder >;
const IndexType rows( 10 ), columns( 10 );
MultidiagonalHost hostMatrix( rows, columns );
- for( IndexType i = 0; i < columns; i++ )
- for( IndexType j = 0; j <= i; j++ )
- hostMatrix.setElement( i, j, i + j );
+ for( IndexType i = 0; i < rows; i++ )
+ for( IndexType j = 0; j < columns; j++ )
+ if( abs( i - j ) <= 1 )
+ hostMatrix.setElement( i, j, i + j );
Matrix matrix( rows, columns );
matrix.getValues() = 0.0;
matrix = hostMatrix;
for( IndexType i = 0; i < columns; i++ )
for( IndexType j = 0; j < rows; j++ )
- {
- if( j > i )
- EXPECT_EQ( matrix.getElement( i, j ), 0.0 );
- else
- EXPECT_EQ( matrix.getElement( i, j ), i + j );
- }
+ if( abs( i - j ) <= 1 )
+ EXPECT_EQ( matrix.getElement( i, j ), i + j );
+ else
+ EXPECT_EQ( matrix.getElement( i, j ), 0.0 );
#ifdef HAVE_CUDA
MultidiagonalCuda cudaMatrix( rows, columns );
- for( IndexType i = 0; i < columns; i++ )
- for( IndexType j = 0; j <= i; j++ )
- cudaMatrix.setElement( i, j, i + j );
+ for( IndexType i = 0; i < rows; i++ )
+ for( IndexType j = 0; j < columns; j++ )
+ if( abs( i - j ) <= 1 )
+ cudaMatrix.setElement( i, j, i + j );
matrix.getValues() = 0.0;
matrix = cudaMatrix;
- for( IndexType i = 0; i < columns; i++ )
- for( IndexType j = 0; j < rows; j++ )
+ for( IndexType i = 0; i < rows; i++ )
+ for( IndexType j = 0; j < columns; j++ )
{
- if( j > i )
- EXPECT_EQ( matrix.getElement( i, j ), 0.0 );
- else
+ if( abs( i - j ) <= 1 )
EXPECT_EQ( matrix.getElement( i, j ), i + j );
+ else
+ EXPECT_EQ( matrix.getElement( i, j ), 0.0 );
}
#endif
}
@@ -1209,123 +1269,125 @@ void test_AssignmentOperator()
template< typename Matrix >
void test_SaveAndLoad()
{
- using RealType = typename Matrix::RealType;
- using DeviceType = typename Matrix::DeviceType;
- using IndexType = typename Matrix::IndexType;
-/*
- * Sets up the following 4x4 dense matrix:
- *
- * / 1 2 3 4 \
- * | 5 6 7 8 |
- * | 9 10 11 12 |
- * \ 13 14 15 16 /
- */
- const IndexType rows = 4;
- const IndexType cols = 4;
+ using RealType = typename Matrix::RealType;
+ using DeviceType = typename Matrix::DeviceType;
+ using IndexType = typename Matrix::IndexType;
+
+ /*
+ * Sets up the following 4x4 matrix:
+ *
+ * / 1 2 0 0 \
+ * | 5 6 7 0 |
+ * | 0 10 11 12 |
+ * \ 0 0 15 16 /
+ */
+ const IndexType rows = 4;
+ const IndexType cols = 4;
- Matrix savedMatrix;
- savedMatrix.reset();
- savedMatrix.setDimensions( rows, cols );
+ Matrix savedMatrix( rows, cols );
- RealType value = 1;
- for( IndexType i = 0; i < rows; i++ )
- for( IndexType j = 0; j < cols; j++ )
- savedMatrix.setElement( i, j, value++ );
-
- ASSERT_NO_THROW( savedMatrix.save( TEST_FILE_NAME ) );
-
- Matrix loadedMatrix;
- loadedMatrix.reset();
- loadedMatrix.setDimensions( rows, cols );
-
- ASSERT_NO_THROW( loadedMatrix.load( TEST_FILE_NAME ) );
-
- EXPECT_EQ( savedMatrix.getElement( 0, 0 ), loadedMatrix.getElement( 0, 0 ) );
- EXPECT_EQ( savedMatrix.getElement( 0, 1 ), loadedMatrix.getElement( 0, 1 ) );
- EXPECT_EQ( savedMatrix.getElement( 0, 2 ), loadedMatrix.getElement( 0, 2 ) );
- EXPECT_EQ( savedMatrix.getElement( 0, 3 ), loadedMatrix.getElement( 0, 3 ) );
-
- EXPECT_EQ( savedMatrix.getElement( 1, 0 ), loadedMatrix.getElement( 1, 0 ) );
- EXPECT_EQ( savedMatrix.getElement( 1, 1 ), loadedMatrix.getElement( 1, 1 ) );
- EXPECT_EQ( savedMatrix.getElement( 1, 2 ), loadedMatrix.getElement( 1, 2 ) );
- EXPECT_EQ( savedMatrix.getElement( 1, 3 ), loadedMatrix.getElement( 1, 3 ) );
-
- EXPECT_EQ( savedMatrix.getElement( 2, 0 ), loadedMatrix.getElement( 2, 0 ) );
- EXPECT_EQ( savedMatrix.getElement( 2, 1 ), loadedMatrix.getElement( 2, 1 ) );
- EXPECT_EQ( savedMatrix.getElement( 2, 2 ), loadedMatrix.getElement( 2, 2 ) );
- EXPECT_EQ( savedMatrix.getElement( 2, 3 ), loadedMatrix.getElement( 2, 3 ) );
-
- EXPECT_EQ( savedMatrix.getElement( 3, 0 ), loadedMatrix.getElement( 3, 0 ) );
- EXPECT_EQ( savedMatrix.getElement( 3, 1 ), loadedMatrix.getElement( 3, 1 ) );
- EXPECT_EQ( savedMatrix.getElement( 3, 2 ), loadedMatrix.getElement( 3, 2 ) );
- EXPECT_EQ( savedMatrix.getElement( 3, 3 ), loadedMatrix.getElement( 3, 3 ) );
-
- EXPECT_EQ( savedMatrix.getElement( 0, 0 ), 1 );
- EXPECT_EQ( savedMatrix.getElement( 0, 1 ), 2 );
- EXPECT_EQ( savedMatrix.getElement( 0, 2 ), 3 );
- EXPECT_EQ( savedMatrix.getElement( 0, 3 ), 4 );
-
- EXPECT_EQ( savedMatrix.getElement( 1, 0 ), 5 );
- EXPECT_EQ( savedMatrix.getElement( 1, 1 ), 6 );
- EXPECT_EQ( savedMatrix.getElement( 1, 2 ), 7 );
- EXPECT_EQ( savedMatrix.getElement( 1, 3 ), 8 );
-
- EXPECT_EQ( savedMatrix.getElement( 2, 0 ), 9 );
- EXPECT_EQ( savedMatrix.getElement( 2, 1 ), 10 );
- EXPECT_EQ( savedMatrix.getElement( 2, 2 ), 11 );
- EXPECT_EQ( savedMatrix.getElement( 2, 3 ), 12 );
-
- EXPECT_EQ( savedMatrix.getElement( 3, 0 ), 13 );
- EXPECT_EQ( savedMatrix.getElement( 3, 1 ), 14 );
- EXPECT_EQ( savedMatrix.getElement( 3, 2 ), 15 );
- EXPECT_EQ( savedMatrix.getElement( 3, 3 ), 16 );
+ RealType value = 1;
+ for( IndexType i = 0; i < rows; i++ )
+ for( IndexType j = 0; j < cols; j++ )
+ {
+ if( abs( i - j ) <= 1 )
+ savedMatrix.setElement( i, j, value );
+ value++;
+ }
+
+ ASSERT_NO_THROW( savedMatrix.save( TEST_FILE_NAME ) );
+
+ Matrix loadedMatrix;
+
+ ASSERT_NO_THROW( loadedMatrix.load( TEST_FILE_NAME ) );
+
+ EXPECT_EQ( savedMatrix.getElement( 0, 0 ), loadedMatrix.getElement( 0, 0 ) );
+ EXPECT_EQ( savedMatrix.getElement( 0, 1 ), loadedMatrix.getElement( 0, 1 ) );
+ EXPECT_EQ( savedMatrix.getElement( 0, 2 ), loadedMatrix.getElement( 0, 2 ) );
+ EXPECT_EQ( savedMatrix.getElement( 0, 3 ), loadedMatrix.getElement( 0, 3 ) );
+
+ EXPECT_EQ( savedMatrix.getElement( 1, 0 ), loadedMatrix.getElement( 1, 0 ) );
+ EXPECT_EQ( savedMatrix.getElement( 1, 1 ), loadedMatrix.getElement( 1, 1 ) );
+ EXPECT_EQ( savedMatrix.getElement( 1, 2 ), loadedMatrix.getElement( 1, 2 ) );
+ EXPECT_EQ( savedMatrix.getElement( 1, 3 ), loadedMatrix.getElement( 1, 3 ) );
+
+ EXPECT_EQ( savedMatrix.getElement( 2, 0 ), loadedMatrix.getElement( 2, 0 ) );
+ EXPECT_EQ( savedMatrix.getElement( 2, 1 ), loadedMatrix.getElement( 2, 1 ) );
+ EXPECT_EQ( savedMatrix.getElement( 2, 2 ), loadedMatrix.getElement( 2, 2 ) );
+ EXPECT_EQ( savedMatrix.getElement( 2, 3 ), loadedMatrix.getElement( 2, 3 ) );
+
+ EXPECT_EQ( savedMatrix.getElement( 3, 0 ), loadedMatrix.getElement( 3, 0 ) );
+ EXPECT_EQ( savedMatrix.getElement( 3, 1 ), loadedMatrix.getElement( 3, 1 ) );
+ EXPECT_EQ( savedMatrix.getElement( 3, 2 ), loadedMatrix.getElement( 3, 2 ) );
+ EXPECT_EQ( savedMatrix.getElement( 3, 3 ), loadedMatrix.getElement( 3, 3 ) );
+
+ EXPECT_EQ( savedMatrix.getElement( 0, 0 ), 1 );
+ EXPECT_EQ( savedMatrix.getElement( 0, 1 ), 2 );
+ EXPECT_EQ( savedMatrix.getElement( 0, 2 ), 0 );
+ EXPECT_EQ( savedMatrix.getElement( 0, 3 ), 0 );
+
+ EXPECT_EQ( savedMatrix.getElement( 1, 0 ), 5 );
+ EXPECT_EQ( savedMatrix.getElement( 1, 1 ), 6 );
+ EXPECT_EQ( savedMatrix.getElement( 1, 2 ), 7 );
+ EXPECT_EQ( savedMatrix.getElement( 1, 3 ), 0 );
+
+ EXPECT_EQ( savedMatrix.getElement( 2, 0 ), 0 );
+ EXPECT_EQ( savedMatrix.getElement( 2, 1 ), 10 );
+ EXPECT_EQ( savedMatrix.getElement( 2, 2 ), 11 );
+ EXPECT_EQ( savedMatrix.getElement( 2, 3 ), 12 );
+
+ EXPECT_EQ( savedMatrix.getElement( 3, 0 ), 0 );
+ EXPECT_EQ( savedMatrix.getElement( 3, 1 ), 0 );
+ EXPECT_EQ( savedMatrix.getElement( 3, 2 ), 15 );
+ EXPECT_EQ( savedMatrix.getElement( 3, 3 ), 16 );
}
template< typename Matrix >
void test_Print()
{
- using RealType = typename Matrix::RealType;
- using DeviceType = typename Matrix::DeviceType;
- using IndexType = typename Matrix::IndexType;
-/*
- * Sets up the following 5x4 sparse matrix:
- *
- * / 1 2 3 4 \
- * | 5 6 7 8 |
- * | 9 10 11 12 |
- * | 13 14 15 16 |
- * \ 17 18 19 20 /
- */
- const IndexType rows = 5;
- const IndexType cols = 4;
+ using RealType = typename Matrix::RealType;
+ using DeviceType = typename Matrix::DeviceType;
+ using IndexType = typename Matrix::IndexType;
- Matrix m;
- m.reset();
- m.setDimensions( rows, cols );
+ /*
+ * Sets up the following 5x4 sparse matrix:
+ *
+ * / 1 2 0 0 \
+ * | 5 6 7 0 |
+ * | 0 10 11 12 |
+ * | 0 0 15 16 |
+ * \ 0 0 0 20 /
+ */
+ const IndexType rows = 5;
+ const IndexType cols = 4;
- RealType value = 1;
- for( IndexType i = 0; i < rows; i++)
- for( IndexType j = 0; j < cols; j++)
- m.setElement( i, j, value++ );
+ Matrix m( rows, cols );
- #include <sstream>
- std::stringstream printed;
- std::stringstream couted;
+ RealType value = 1;
+ for( IndexType i = 0; i < rows; i++)
+ for( IndexType j = 0; j < cols; j++)
+ {
+ if( abs( i - j ) <= 1 )
+ m.setElement( i, j, value );
+ value++;
+ }
- //change the underlying buffer and save the old buffer
- auto old_buf = std::cout.rdbuf(printed.rdbuf());
+ std::stringstream printed;
+ std::stringstream couted;
- m.print( std::cout ); //all the std::cout goes to ss
+ //change the underlying buffer and save the old buffer
+ auto old_buf = std::cout.rdbuf(printed.rdbuf());
- std::cout.rdbuf(old_buf); //reset
+ m.print( std::cout ); //all the std::cout goes to ss
- couted << "Row: 0 -> Col:0->1 Col:1->2 Col:2->3 Col:3->4\t\n"
- "Row: 1 -> Col:0->5 Col:1->6 Col:2->7 Col:3->8\t\n"
- "Row: 2 -> Col:0->9 Col:1->10 Col:2->11 Col:3->12\t\n"
- "Row: 3 -> Col:0->13 Col:1->14 Col:2->15 Col:3->16\t\n"
- "Row: 4 -> Col:0->17 Col:1->18 Col:2->19 Col:3->20\t\n";
+ std::cout.rdbuf(old_buf); //reset
+ couted << "Row: 0 -> Col:0->1\t Col:1->2\t\n"
+ "Row: 1 -> Col:0->5\t Col:1->6\t Col:2->7\t\n"
+ "Row: 2 -> Col:1->10\t Col:2->11\t Col:3->12\t\n"
+ "Row: 3 -> Col:2->15\t Col:3->16\t\n"
+ "Row: 4 -> Col:3->20\t\n";
- EXPECT_EQ( printed.str(), couted.str() );
+ EXPECT_EQ( printed.str(), couted.str() );
}
// test fixture for typed tests
@@ -1388,6 +1450,21 @@ TYPED_TEST( MatrixTest, setLikeTest )
test_SetLike< MatrixType, MatrixType >();
}
+TYPED_TEST( MatrixTest, getNonemptyRowsCountTest )
+{
+ using MatrixType = typename TestFixture::MatrixType;
+
+ test_GetNonemptyRowsCount< MatrixType >();
+}
+
+/*
+TYPED_TEST( MatrixTest, getCompressedRowLengthTest )
+{
+ using MatrixType = typename TestFixture::MatrixType;
+
+ test_GetCompressedRowLengths< MatrixType >();
+}
+
TYPED_TEST( MatrixTest, getRowLengthTest )
{
using MatrixType = typename TestFixture::MatrixType;
@@ -1395,11 +1472,11 @@ TYPED_TEST( MatrixTest, getRowLengthTest )
test_GetRowLength< MatrixType >();
}
-TYPED_TEST( MatrixTest, getNumberOfMatrixElementsTest )
+TYPED_TEST( MatrixTest, getAllocatedElementsCountTest )
{
using MatrixType = typename TestFixture::MatrixType;
- test_GetNumberOfMatrixElements< MatrixType >();
+ test_GetAllocatedElementsCount< MatrixType >();
}
TYPED_TEST( MatrixTest, getNumberOfNonzeroMatrixElementsTest )
@@ -1465,6 +1542,19 @@ TYPED_TEST( MatrixTest, addMatrixTest )
test_AddMatrix< MatrixType >();
}
+TYPED_TEST( MatrixTest, addMatrixTest_differentOrdering )
+{
+ using MatrixType = typename TestFixture::MatrixType;
+
+ using RealType = typename MatrixType::RealType;
+ using DeviceType = typename MatrixType::DeviceType;
+ using IndexType = typename MatrixType::IndexType;
+ using RealAllocatorType = typename MatrixType::RealAllocatorType;
+ using MatrixType2 = TNL::Matrices::Multidiagonal< RealType, DeviceType, IndexType, ! MatrixType::getRowMajorOrder(), RealAllocatorType >;
+
+ test_AddMatrix< MatrixType, MatrixType2 >();
+}
+
TYPED_TEST( MatrixTest, assignmentOperatorTest )
{
using MatrixType = typename TestFixture::MatrixType;
@@ -1485,6 +1575,7 @@ TYPED_TEST( MatrixTest, printTest )
test_Print< MatrixType >();
}
+*/
//// test_getType is not general enough yet. DO NOT TEST IT YET.