From b20418121b70638c87ce82851777105c29e9e464 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= <oberhuber.tomas@gmail.com>
Date: Mon, 13 Jan 2020 21:39:52 +0100
Subject: [PATCH] Debugging multidiagonal matrix.
---
src/TNL/Matrices/Matrix.h | 4 +-
src/TNL/Matrices/MatrixView.h | 4 +-
src/TNL/Matrices/Multidiagonal.h | 10 +-
src/TNL/Matrices/Multidiagonal.hpp | 19 +-
src/TNL/Matrices/MultidiagonalMatrixRowView.h | 11 +-
.../Matrices/MultidiagonalMatrixRowView.hpp | 35 +-
src/TNL/Matrices/MultidiagonalMatrixView.h | 7 +-
src/TNL/Matrices/MultidiagonalMatrixView.hpp | 172 +++--
src/TNL/Matrices/Tridiagonal.h | 4 +-
src/TNL/Matrices/Tridiagonal.hpp | 17 +-
src/TNL/Matrices/TridiagonalMatrixView.h | 4 +-
src/TNL/Matrices/TridiagonalMatrixView.hpp | 10 +-
.../details/MultidiagonalMatrixIndexer.h | 3 +
.../Matrices/MultidiagonalMatrixTest.h | 666 ++++++++----------
14 files changed, 453 insertions(+), 513 deletions(-)
diff --git a/src/TNL/Matrices/Matrix.h b/src/TNL/Matrices/Matrix.h
index 7813fa9625..ebe7ccc21f 100644
--- a/src/TNL/Matrices/Matrix.h
+++ b/src/TNL/Matrices/Matrix.h
@@ -76,11 +76,11 @@ public:
__cuda_callable__
IndexType getColumns() const;
- virtual bool setElement( const IndexType row,
+ virtual void setElement( const IndexType row,
const IndexType column,
const RealType& value ) = 0;
- virtual bool addElement( const IndexType row,
+ virtual void addElement( const IndexType row,
const IndexType column,
const RealType& value,
const RealType& thisElementMultiplicator = 1.0 ) = 0;
diff --git a/src/TNL/Matrices/MatrixView.h b/src/TNL/Matrices/MatrixView.h
index 467d023496..2a6429df58 100644
--- a/src/TNL/Matrices/MatrixView.h
+++ b/src/TNL/Matrices/MatrixView.h
@@ -73,11 +73,11 @@ public:
* in the future and it does not slow down, declare them as virtual here.
*/
- virtual bool setElement( const IndexType row,
+ virtual void setElement( const IndexType row,
const IndexType column,
const RealType& value ) = 0;
- virtual bool addElement( const IndexType row,
+ virtual void addElement( const IndexType row,
const IndexType column,
const RealType& value,
const RealType& thisElementMultiplicator = 1.0 ) = 0;
diff --git a/src/TNL/Matrices/Multidiagonal.h b/src/TNL/Matrices/Multidiagonal.h
index 5d23cd960e..1741c0c75f 100644
--- a/src/TNL/Matrices/Multidiagonal.h
+++ b/src/TNL/Matrices/Multidiagonal.h
@@ -38,12 +38,12 @@ class Multidiagonal : public 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 DiagonalsShiftsType = Containers::Vector< IndexType, DeviceType, IndexType, IndexAllocatorType >;
+ using DiagonalsShiftsView = typename DiagonalsShiftsType::ViewType;
+ using RowView = MultidiagonalMatrixRowView< ValuesViewType, IndexerType, DiagonalsShiftsView >;
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;
@@ -119,11 +119,11 @@ class Multidiagonal : public Matrix< Real, Device, Index, RealAllocator >
void setValue( const RealType& v );
- bool setElement( const IndexType row,
+ void setElement( const IndexType row,
const IndexType column,
const RealType& value );
- bool addElement( const IndexType row,
+ void addElement( const IndexType row,
const IndexType column,
const RealType& value,
const RealType& thisElementMultiplicator = 1.0 );
diff --git a/src/TNL/Matrices/Multidiagonal.hpp b/src/TNL/Matrices/Multidiagonal.hpp
index 95f6667c1c..53e3c7f2f4 100644
--- a/src/TNL/Matrices/Multidiagonal.hpp
+++ b/src/TNL/Matrices/Multidiagonal.hpp
@@ -44,6 +44,7 @@ Multidiagonal( const IndexType rows,
const IndexType columns,
const Vector& diagonalsShifts )
{
+ TNL_ASSERT_GT( diagonalsShifts.getSize(), 0, "Cannot construct mutltidiagonal matrix with no diagonals shifts." );
this->setDimensions( rows, columns, diagonalsShifts );
}
@@ -60,6 +61,7 @@ getView() const -> ViewType
// TODO: fix when getConstView works
return ViewType( const_cast< Multidiagonal* >( this )->values.getView(),
const_cast< Multidiagonal* >( this )->diagonalsShifts.getView(),
+ const_cast< Multidiagonal* >( this )->hostDiagonalsShifts.getView(),
indexer );
}
@@ -358,11 +360,11 @@ template< typename Real,
bool RowMajorOrder,
typename RealAllocator,
typename IndexAllocator >
-bool
+void
Multidiagonal< Real, Device, Index, RowMajorOrder, RealAllocator, IndexAllocator >::
setElement( const IndexType row, const IndexType column, const RealType& value )
{
- return this->view.setElement( row, column, value );
+ this->view.setElement( row, column, value );
}
template< typename Real,
@@ -371,14 +373,14 @@ template< typename Real,
bool RowMajorOrder,
typename RealAllocator,
typename IndexAllocator >
-bool
+void
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 );
+ this->view.addElement( row, column, value, thisElementMultiplicator );
}
template< typename Real,
@@ -745,14 +747,7 @@ 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;
- }
+ this->view.print( str );
}
template< typename Real,
diff --git a/src/TNL/Matrices/MultidiagonalMatrixRowView.h b/src/TNL/Matrices/MultidiagonalMatrixRowView.h
index 68b5be55c1..0825d6fb36 100644
--- a/src/TNL/Matrices/MultidiagonalMatrixRowView.h
+++ b/src/TNL/Matrices/MultidiagonalMatrixRowView.h
@@ -14,7 +14,8 @@ namespace TNL {
namespace Matrices {
template< typename ValuesView,
- typename Indexer >
+ typename Indexer,
+ typename DiagonalsShiftsView_ >
class MultidiagonalMatrixRowView
{
public:
@@ -23,11 +24,13 @@ class MultidiagonalMatrixRowView
using IndexType = typename ValuesView::IndexType;
using ValuesViewType = ValuesView;
using IndexerType = Indexer;
+ using DiagonalsShiftsView = DiagonalsShiftsView_;
__cuda_callable__
MultidiagonalMatrixRowView( const IndexType rowIdx,
- const ValuesViewType& values,
- const IndexerType& indexer );
+ const DiagonalsShiftsView& diagonalsShifts,
+ const ValuesViewType& values,
+ const IndexerType& indexer);
__cuda_callable__
IndexType getSize() const;
@@ -48,6 +51,8 @@ class MultidiagonalMatrixRowView
IndexType rowIdx;
+ DiagonalsShiftsView diagonalsShifts;
+
ValuesViewType values;
Indexer indexer;
diff --git a/src/TNL/Matrices/MultidiagonalMatrixRowView.hpp b/src/TNL/Matrices/MultidiagonalMatrixRowView.hpp
index 349fbe8ea3..88aae3f150 100644
--- a/src/TNL/Matrices/MultidiagonalMatrixRowView.hpp
+++ b/src/TNL/Matrices/MultidiagonalMatrixRowView.hpp
@@ -13,58 +13,59 @@
namespace TNL {
namespace Matrices {
-template< typename ValuesView, typename Indexer >
+template< typename ValuesView, typename Indexer, typename DiagonalsShiftsView >
__cuda_callable__
-MultidiagonalMatrixRowView< ValuesView, Indexer >::
+MultidiagonalMatrixRowView< ValuesView, Indexer, DiagonalsShiftsView >::
MultidiagonalMatrixRowView( const IndexType rowIdx,
- const ValuesViewType& values,
- const IndexerType& indexer )
-: rowIdx( rowIdx ), values( values ), indexer( indexer )
+ const DiagonalsShiftsView& diagonalsShifts,
+ const ValuesViewType& values,
+ const IndexerType& indexer )
+: rowIdx( rowIdx ), diagonalsShifts( diagonalsShifts ), values( values ), indexer( indexer )
{
}
-template< typename ValuesView, typename Indexer >
+template< typename ValuesView, typename Indexer, typename DiagonalsShiftsView >
__cuda_callable__
auto
-MultidiagonalMatrixRowView< ValuesView, Indexer >::
+MultidiagonalMatrixRowView< ValuesView, Indexer, DiagonalsShiftsView >::
getSize() const -> IndexType
{
return indexer.getRowSize();
}
-template< typename ValuesView, typename Indexer >
+template< typename ValuesView, typename Indexer, typename DiagonalsShiftsView >
__cuda_callable__
auto
-MultidiagonalMatrixRowView< ValuesView, Indexer >::
+MultidiagonalMatrixRowView< ValuesView, Indexer, DiagonalsShiftsView >::
getColumnIndex( const IndexType localIdx ) const -> const IndexType
{
TNL_ASSERT_GE( localIdx, 0, "" );
- TNL_ASSERT_LT( localIdx, 3, "" );
- return rowIdx + localIdx - 1;
+ TNL_ASSERT_LT( localIdx, indexer.getDiagonals(), "" );
+ return rowIdx + diagonalsShifts[ localIdx ];
}
-template< typename ValuesView, typename Indexer >
+template< typename ValuesView, typename Indexer, typename DiagonalsShiftsView >
__cuda_callable__
auto
-MultidiagonalMatrixRowView< ValuesView, Indexer >::
+MultidiagonalMatrixRowView< ValuesView, Indexer, DiagonalsShiftsView >::
getValue( const IndexType localIdx ) const -> const RealType&
{
return this->values[ this->indexer.getGlobalIndex( rowIdx, localIdx ) ];
}
-template< typename ValuesView, typename Indexer >
+template< typename ValuesView, typename Indexer, typename DiagonalsShiftsView >
__cuda_callable__
auto
-MultidiagonalMatrixRowView< ValuesView, Indexer >::
+MultidiagonalMatrixRowView< ValuesView, Indexer, DiagonalsShiftsView >::
getValue( const IndexType localIdx ) -> RealType&
{
return this->values[ this->indexer.getGlobalIndex( rowIdx, localIdx ) ];
}
-template< typename ValuesView, typename Indexer >
+template< typename ValuesView, typename Indexer, typename DiagonalsShiftsView >
__cuda_callable__
void
-MultidiagonalMatrixRowView< ValuesView, Indexer >::
+MultidiagonalMatrixRowView< ValuesView, Indexer, DiagonalsShiftsView >::
setElement( const IndexType localIdx,
const RealType& value )
{
diff --git a/src/TNL/Matrices/MultidiagonalMatrixView.h b/src/TNL/Matrices/MultidiagonalMatrixView.h
index addeb18b34..3d33ac0aea 100644
--- a/src/TNL/Matrices/MultidiagonalMatrixView.h
+++ b/src/TNL/Matrices/MultidiagonalMatrixView.h
@@ -38,7 +38,7 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
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 >;
+ using RowView = MultidiagonalMatrixRowView< ValuesViewType, IndexerType, DiagonalsShiftsView >;
// TODO: remove this - it is here only for compatibility with original matrix implementation
typedef Containers::Vector< IndexType, DeviceType, IndexType > CompressedRowLengthsVector;
@@ -55,6 +55,7 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
MultidiagonalMatrixView( const ValuesViewType& values,
const DiagonalsShiftsView& diagonalsShifts,
+ const HostDiagonalsShiftsView& hostDiagonalsShifts,
const IndexerType& indexer );
ViewType getView();
@@ -92,11 +93,11 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
void setValue( const RealType& v );
- bool setElement( const IndexType row,
+ void setElement( const IndexType row,
const IndexType column,
const RealType& value );
- bool addElement( const IndexType row,
+ void addElement( const IndexType row,
const IndexType column,
const RealType& value,
const RealType& thisElementMultiplicator = 1.0 );
diff --git a/src/TNL/Matrices/MultidiagonalMatrixView.hpp b/src/TNL/Matrices/MultidiagonalMatrixView.hpp
index 3d9b0237f9..1ba8dc34d0 100644
--- a/src/TNL/Matrices/MultidiagonalMatrixView.hpp
+++ b/src/TNL/Matrices/MultidiagonalMatrixView.hpp
@@ -33,9 +33,11 @@ template< typename Real,
MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
MultidiagonalMatrixView( const ValuesViewType& values,
const DiagonalsShiftsView& diagonalsShifts,
+ const HostDiagonalsShiftsView& hostDiagonalsShifts,
const IndexerType& indexer )
: MatrixView< Real, Device, Index >( indexer.getRows(), indexer.getColumns(), values ),
diagonalsShifts( diagonalsShifts ),
+ hostDiagonalsShifts( hostDiagonalsShifts ),
indexer( indexer )
{
}
@@ -48,7 +50,10 @@ auto
MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
getView() -> ViewType
{
- return ViewType( this->values.getView(), indexer );
+ return ViewType( const_cast< MultidiagonalMatrixView* >( this )->values.getView(),
+ const_cast< MultidiagonalMatrixView* >( this )->diagonalsShifts.getView(),
+ const_cast< MultidiagonalMatrixView* >( this )->hostDiagonalsShifts.getView(),
+ indexer );
}
template< typename Real,
@@ -59,7 +64,10 @@ auto
MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
getConstView() const -> ConstViewType
{
- return ConstViewType( this->values.getConstView(), indexer );
+ return ConstViewType( this->values.getConstView(),
+ this->diagonalsShifts.getConstView(),
+ this->hostDiagonalsShifts.getConstView(),
+ indexer );
}
template< typename Real,
@@ -208,7 +216,11 @@ void
MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
setValue( const RealType& v )
{
- this->values = v;
+ const RealType newValue = v;
+ auto f = [=] __cuda_callable__ ( const IndexType& rowIdx, const IndexType& localIdx, const IndexType columnIdx, RealType& value ) mutable {
+ value = newValue;
+ };
+ this->forAllRows( f );
}
template< typename Real,
@@ -220,7 +232,7 @@ auto
MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
getRow( const IndexType& rowIdx ) const -> const RowView
{
- return RowView( rowIdx, this->values.getView(), this->indexer );
+ return RowView( rowIdx, this->diagonalsShifts.getView(), this->values.getView(), this->indexer );
}
template< typename Real,
@@ -232,14 +244,14 @@ auto
MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
getRow( const IndexType& rowIdx ) -> RowView
{
- return RowView( rowIdx, this->values.getView(), this->indexer );
+ return RowView( rowIdx, this->diagonalsShifts.getView(), this->values.getView(), this->indexer );
}
template< typename Real,
typename Device,
typename Index,
bool RowMajorOrder >
-bool
+void
MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
setElement( const IndexType row, const IndexType column, const RealType& value )
{
@@ -247,21 +259,26 @@ setElement( const IndexType row, const IndexType column, const RealType& value )
TNL_ASSERT_LT( row, this->getRows(), "" );
TNL_ASSERT_GE( column, 0, "" );
TNL_ASSERT_LT( column, this->getColumns(), "" );
- if( abs( row - column ) > 1 )
+
+ for( IndexType i = 0; i < hostDiagonalsShifts.getSize(); i++ )
+ if( row + hostDiagonalsShifts[ i ] == column )
+ {
+ this->values.setElement( this->getElementIndex( row, i ), value );
+ return;
+ }
+ if( value != 0.0 )
{
std::stringstream msg;
- msg << "Wrong matrix element coordinates ( " << row << ", " << column << " ) in tridiagonal matrix.";
+ msg << "Wrong matrix element coordinates ( " << row << ", " << column << " ) in multidiagonal 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
+void
MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
addElement( const IndexType row,
const IndexType column,
@@ -272,15 +289,20 @@ addElement( const IndexType row,
TNL_ASSERT_LT( row, this->getRows(), "" );
TNL_ASSERT_GE( column, 0, "" );
TNL_ASSERT_LT( column, this->getColumns(), "" );
- if( abs( row - column ) > 1 )
+
+ for( IndexType i = 0; i < hostDiagonalsShifts.getSize(); i++ )
+ if( row + hostDiagonalsShifts[ i ] == column )
+ {
+ const Index idx = this->getElementIndex( row, i );
+ this->values.setElement( idx, thisElementMultiplicator * this->values.getElement( idx ) + value );
+ return;
+ }
+ if( value != 0.0 )
{
std::stringstream msg;
- msg << "Wrong matrix element coordinates ( " << row << ", " << column << " ) in tridiagonal matrix.";
+ msg << "Wrong matrix element coordinates ( " << row << ", " << column << " ) in multidiagonal 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,
@@ -296,9 +318,10 @@ getElement( const IndexType row, const IndexType column ) const
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 ) );
+ for( IndexType i = 0; i < hostDiagonalsShifts.getSize(); i++ )
+ if( row + hostDiagonalsShifts[ i ] == column )
+ return this->values.getElement( this->getElementIndex( row, i ) );
+ return 0.0;
}
template< typename Real,
@@ -326,35 +349,20 @@ rowsReduction( IndexType first, IndexType last, Fetch& fetch, Reduce& reduce, Ke
{
using Real_ = decltype( fetch( IndexType(), IndexType(), RealType() ) );
const auto values_view = this->values.getConstView();
+ const auto diagonalsShifts_view = this->diagonalsShifts.getConstView();
+ const IndexType diagonalsCount = this->diagonalsShifts.getSize();
+ const IndexType columns = this->getColumns();
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() )
+ for( IndexType localIdx = 0; localIdx < diagonalsCount; localIdx++ )
{
- 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 ) ] ) );
+ const IndexType columnIdx = rowIdx + diagonalsShifts_view[ localIdx ];
+ if( columnIdx >= 0 && columnIdx < columns )
+ reduce( sum, fetch( rowIdx, columnIdx, values_view[ indexer.getGlobalIndex( rowIdx, localIdx ) ] ) );
}
+ keep( rowIdx, sum );
};
Algorithms::ParallelFor< DeviceType >::exec( first, last, f );
}
@@ -368,7 +376,7 @@ 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 );
+ this->rowsReduction( 0, this->indexer.getNonemptyRowsCount(), fetch, reduce, keep, zero );
}
template< typename Real,
@@ -381,26 +389,17 @@ MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
forRows( IndexType first, IndexType last, Function& function ) const
{
const auto values_view = this->values.getConstView();
+ const auto diagonalsShifts_view = this->diagonalsShifts.getConstView();
+ const IndexType diagonalsCount = this->diagonalsShifts.getSize();
+ const IndexType columns = this->getColumns();
const auto indexer = this->indexer;
auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
- if( rowIdx == 0 )
+ for( IndexType localIdx = 0; localIdx < diagonalsCount; localIdx++ )
{
- 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 ) ] );
+ const IndexType columnIdx = rowIdx + diagonalsShifts_view[ localIdx ];
+ if( columnIdx >= 0 && columnIdx < columns )
+ function( rowIdx, localIdx, columnIdx, values_view[ indexer.getGlobalIndex( rowIdx, localIdx, 0 ) ] );
}
- 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 );
}
@@ -415,26 +414,17 @@ MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
forRows( IndexType first, IndexType last, Function& function )
{
auto values_view = this->values.getView();
+ const auto diagonalsShifts_view = this->diagonalsShifts.getConstView();
+ const IndexType diagonalsCount = this->diagonalsShifts.getSize();
+ const IndexType columns = this->getColumns();
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() )
+ for( IndexType localIdx = 0; localIdx < diagonalsCount; localIdx++ )
{
- function( rowIdx, 0, rowIdx - 1, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ] );
- function( rowIdx, 1, rowIdx, values_view[ indexer.getGlobalIndex( rowIdx, 1 ) ] );
+ const IndexType columnIdx = rowIdx + diagonalsShifts_view[ localIdx ];
+ if( columnIdx >= 0 && columnIdx < columns )
+ function( rowIdx, localIdx, columnIdx, values_view[ indexer.getGlobalIndex( rowIdx, localIdx ) ] );
}
- else
- function( rowIdx, 0, rowIdx, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ] );
};
Algorithms::ParallelFor< DeviceType >::exec( first, last, f );
}
@@ -460,7 +450,7 @@ void
MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
forAllRows( Function& function )
{
- this->forRows( 0, this->indexer.getNonEmptyRowsCount(), function );
+ this->forRows( 0, this->indexer.getNonemptyRowsCount(), function );
}
template< typename Real,
@@ -517,7 +507,7 @@ addMatrix( const MultidiagonalMatrixView< Real_, Device_, Index_, RowMajorOrder_
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( RowMajorOrder == RowMajorOrder_ )
{
if( thisMatrixMultiplicator == 1.0 )
this->values += matrixMultiplicator * matrix.getValues();
@@ -544,7 +534,7 @@ addMatrix( const MultidiagonalMatrixView< Real_, Device_, Index_, RowMajorOrder_
this->forAllRows( add1 );
else
this->forAllRows( addGen );
- }
+ }*/
}
#ifdef HAVE_CUDA
@@ -673,12 +663,19 @@ template< typename Real,
bool RowMajorOrder >
void MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::print( std::ostream& str ) const
{
- for( IndexType row = 0; row < this->getRows(); row++ )
+ for( IndexType rowIdx = 0; rowIdx < this->getRows(); rowIdx++ )
{
- 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 <<"Row: " << rowIdx << " -> ";
+ for( IndexType localIdx = 0; localIdx < this->hostDiagonalsShifts.getSize(); localIdx++ )
+ {
+ const IndexType columnIdx = rowIdx + this->hostDiagonalsShifts[ localIdx ];
+ if( columnIdx >= 0 && columnIdx < this->columns )
+ {
+ auto v = this->values.getElement( this->indexer.getGlobalIndex( rowIdx, localIdx ) );
+ if( v )
+ str << " Col:" << columnIdx << "->" << v << "\t";
+ }
+ }
str << std::endl;
}
}
@@ -713,15 +710,8 @@ template< typename Real,
bool RowMajorOrder >
__cuda_callable__
Index MultidiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
-getElementIndex( const IndexType row, const IndexType column ) const
+getElementIndex( const IndexType row, const IndexType localIdx ) 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 );
}
diff --git a/src/TNL/Matrices/Tridiagonal.h b/src/TNL/Matrices/Tridiagonal.h
index e7e3ab6b27..82549e7445 100644
--- a/src/TNL/Matrices/Tridiagonal.h
+++ b/src/TNL/Matrices/Tridiagonal.h
@@ -97,11 +97,11 @@ class Tridiagonal : public Matrix< Real, Device, Index, RealAllocator >
void setValue( const RealType& v );
- bool setElement( const IndexType row,
+ void setElement( const IndexType row,
const IndexType column,
const RealType& value );
- bool addElement( const IndexType row,
+ void addElement( const IndexType row,
const IndexType column,
const RealType& value,
const RealType& thisElementMultiplicator = 1.0 );
diff --git a/src/TNL/Matrices/Tridiagonal.hpp b/src/TNL/Matrices/Tridiagonal.hpp
index 6c09238ff0..41d722c6a7 100644
--- a/src/TNL/Matrices/Tridiagonal.hpp
+++ b/src/TNL/Matrices/Tridiagonal.hpp
@@ -285,11 +285,11 @@ template< typename Real,
typename Index,
bool RowMajorOrder,
typename RealAllocator >
-bool
+void
Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
setElement( const IndexType row, const IndexType column, const RealType& value )
{
- return this->view.setElement( row, column, value );
+ this->view.setElement( row, column, value );
}
template< typename Real,
@@ -297,14 +297,14 @@ template< typename Real,
typename Index,
bool RowMajorOrder,
typename RealAllocator >
-bool
+void
Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
addElement( const IndexType row,
const IndexType column,
const RealType& value,
const RealType& thisElementMultiplicator )
{
- return this->view.addElement( row, column, value, thisElementMultiplicator );
+ this->view.addElement( row, column, value, thisElementMultiplicator );
}
template< typename Real,
@@ -645,14 +645,7 @@ void
Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
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;
- }
+ this->view.print( str );
}
template< typename Real,
diff --git a/src/TNL/Matrices/TridiagonalMatrixView.h b/src/TNL/Matrices/TridiagonalMatrixView.h
index 128b484947..7db517cbd3 100644
--- a/src/TNL/Matrices/TridiagonalMatrixView.h
+++ b/src/TNL/Matrices/TridiagonalMatrixView.h
@@ -81,11 +81,11 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
void setValue( const RealType& v );
- bool setElement( const IndexType row,
+ void setElement( const IndexType row,
const IndexType column,
const RealType& value );
- bool addElement( const IndexType row,
+ void addElement( const IndexType row,
const IndexType column,
const RealType& value,
const RealType& thisElementMultiplicator = 1.0 );
diff --git a/src/TNL/Matrices/TridiagonalMatrixView.hpp b/src/TNL/Matrices/TridiagonalMatrixView.hpp
index 4d4950c4e9..e851d2a1f9 100644
--- a/src/TNL/Matrices/TridiagonalMatrixView.hpp
+++ b/src/TNL/Matrices/TridiagonalMatrixView.hpp
@@ -213,7 +213,7 @@ template< typename Real,
typename Device,
typename Index,
bool RowMajorOrder >
-bool
+void
TridiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
setElement( const IndexType row, const IndexType column, const RealType& value )
{
@@ -235,7 +235,7 @@ template< typename Real,
typename Device,
typename Index,
bool RowMajorOrder >
-bool
+void
TridiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
addElement( const IndexType row,
const IndexType column,
@@ -638,7 +638,11 @@ void TridiagonalMatrixView< Real, Device, Index, RowMajorOrder >::print( std::os
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";
+ {
+ auto v = this->getElement( row, column );
+ if( v )
+ str << " Col:" << column << "->" << v << "\t";
+ }
str << std::endl;
}
}
diff --git a/src/TNL/Matrices/details/MultidiagonalMatrixIndexer.h b/src/TNL/Matrices/details/MultidiagonalMatrixIndexer.h
index 0f0436d748..3597c30f7d 100644
--- a/src/TNL/Matrices/details/MultidiagonalMatrixIndexer.h
+++ b/src/TNL/Matrices/details/MultidiagonalMatrixIndexer.h
@@ -77,6 +77,9 @@ class MultidiagonalMatrixIndexer
__cuda_callable__
const IndexType& getColumns() const { return this->columns; };
+ __cuda_callable__
+ const IndexType& getDiagonals() const { return this->diagonals; };
+
__cuda_callable__
const IndexType& getNonemptyRowsCount() const { return this->nonemptyRows; };
diff --git a/src/UnitTests/Matrices/MultidiagonalMatrixTest.h b/src/UnitTests/Matrices/MultidiagonalMatrixTest.h
index cb9916e4c9..514ea39e0e 100644
--- a/src/UnitTests/Matrices/MultidiagonalMatrixTest.h
+++ b/src/UnitTests/Matrices/MultidiagonalMatrixTest.h
@@ -137,7 +137,6 @@ void test_GetNonemptyRowsCount()
Matrix m3( 8, 5, DiagonalsShiftsType({ -2, 0, 3, 5 }) );
m3.setValue( 1.0 );
EXPECT_EQ( m3.getNonemptyRowsCount(), 7 );
-
}
template< typename Matrix >
@@ -148,87 +147,53 @@ void test_GetCompressedRowLengths()
using IndexType = typename Matrix::IndexType;
using DiagonalsShiftsType = typename Matrix::DiagonalsShiftsType;
- const IndexType rows = 10;
- const IndexType cols = 11;
-
- Matrix m( rows, cols );
-
- // Insert values into the rows.
- RealType value = 1;
-
- 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 -> 3 elements
- m.setElement( 1, i, value++ );
-
- for( IndexType i = 1; i < 3; i++ ) // 2nd row -> 2 elements
- m.setElement( 2, i, value++ );
-
- for( IndexType i = 2; i < 5; i++ ) // 3rd row -> 3 elements
- m.setElement( 3, i, value++ );
-
- for( IndexType i = 3; i < 6; i++ ) // 4th row -> 3 elements
- m.setElement( 4, i, value++ );
-
- for( IndexType i = 4; i < 6; i++ ) // 5th row -> 2 elements
- m.setElement( 5, i, value++ );
-
- for( IndexType i = 5; i < 8; i++ ) // 6th row -> 3 elements
- m.setElement( 6, i, value++ );
-
- for( IndexType i = 6; i < 8; i++ ) // 7th row -> 2 elements
- m.setElement( 7, i, value++ );
-
- for( IndexType i = 7; i < 10; i++ ) // 8th row -> 3 elements
- m.setElement( 8, i, value++ );
+ /*
+ * Sets up the following 8x8 matrix:
+ *
+ * / 0 0 0 1 0 1 0 0 \ -> 2
+ * | 0 1 0 0 1 0 1 0 | -> 3
+ * | 1 0 1 0 0 1 0 1 | -> 4
+ * | 0 1 0 1 0 0 1 0 | -> 3
+ * | 0 0 1 0 1 0 0 1 | -> 3
+ * | 0 0 0 1 0 1 0 0 | -> 2
+ * | 0 0 0 0 1 0 1 0 | -> 2
+ * \ 0 0 0 0 0 1 0 0 / -> 1
+ */
+
+ const IndexType rows = 8;
+ const IndexType cols = 8;
- for( IndexType i = 8; i < 11; i++ ) // 9th row -> 3 elements
- m.setElement( 9, i, value++ );
+ Matrix m( rows, cols, DiagonalsShiftsType({ -2, 0, 3, 5 }) );
+ m.setValue( 1.0 );
+ m.setElement( 0, 0, 0.0 );
+ m.setElement( 7, 7, 0.0 );
typename Matrix::CompressedRowLengthsVector rowLengths( rows );
rowLengths = 0;
m.getCompressedRowLengths( rowLengths );
- typename Matrix::CompressedRowLengthsVector correctRowLengths{ 2, 3, 2, 3, 3, 2, 3, 2, 3, 3 };
+ typename Matrix::CompressedRowLengthsVector correctRowLengths{ 2, 3, 4, 3, 3, 2, 2, 1 };
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 DiagonalsShiftsType = typename Matrix::DiagonalsShiftsType;
-
- const IndexType rows = 8;
- const IndexType cols = 7;
-
- Matrix m( rows, cols );
-
- 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_GetAllocatedElementsCount()
{
using RealType = typename Matrix::RealType;
using DeviceType = typename Matrix::DeviceType;
using IndexType = typename Matrix::IndexType;
+ using DiagonalsShiftsType = typename Matrix::DiagonalsShiftsType;
const IndexType rows = 7;
const IndexType cols = 6;
- Matrix m( rows, cols );
+ Matrix m1( 7, 6, DiagonalsShiftsType( { -2, 0, 3, 5 } ) );
+ EXPECT_EQ( m1.getAllocatedElementsCount(), 28 );
+
+ Matrix m2( 8, 6, DiagonalsShiftsType( { -2, 0, 3, 5 } ) );
+ EXPECT_EQ( m2.getAllocatedElementsCount(), 32 );
- EXPECT_EQ( m.getAllocatedElementsCount(), 21 );
+ Matrix m3( 9, 6, DiagonalsShiftsType( { -2, 0, 3, 5 } ) );
+ EXPECT_EQ( m3.getAllocatedElementsCount(), 32 );
}
template< typename Matrix >
@@ -237,29 +202,27 @@ void test_GetNumberOfNonzeroMatrixElements()
using RealType = typename Matrix::RealType;
using DeviceType = typename Matrix::DeviceType;
using IndexType = typename Matrix::IndexType;
+ using DiagonalsShiftsType = typename Matrix::DiagonalsShiftsType;
/*
* 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 /
+ * / 0 0 1 0 1 0 \ -> 2
+ * | 0 1 0 1 0 1 | -> 3
+ * | 0 0 1 0 1 0 | -> 2
+ * | 1 0 0 1 0 1 | -> 3
+ * | 0 1 0 0 1 0 | -> 2
+ * | 0 0 1 0 0 1 | -> 2
+ * \ 0 0 0 1 0 0 / -> 1
+ * ----
+ * 15
*/
const IndexType rows = 7;
const IndexType cols = 6;
- Matrix m( rows, cols );
-
- 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);
+ Matrix m( rows, cols, DiagonalsShiftsType( { -3, 0, 2, 4 } ) );
+ m.setValue( 1.0 );
+ m.setElement( 0, 0, 0.0 );
EXPECT_EQ( m.getNumberOfNonzeroMatrixElements(), 15 );
}
@@ -270,6 +233,7 @@ void test_Reset()
using RealType = typename Matrix::RealType;
using DeviceType = typename Matrix::DeviceType;
using IndexType = typename Matrix::IndexType;
+ using DiagonalsShiftsType = typename Matrix::DiagonalsShiftsType;
/*
* Sets up the following 5x4 matrix:
@@ -283,8 +247,7 @@ void test_Reset()
const IndexType rows = 5;
const IndexType cols = 4;
- Matrix m( rows, cols );
-
+ Matrix m( rows, cols, DiagonalsShiftsType( { 0, 1, 2, 4 } ) );
m.reset();
EXPECT_EQ( m.getRows(), 0 );
@@ -297,130 +260,73 @@ void test_SetValue()
using RealType = typename Matrix::RealType;
using DeviceType = typename Matrix::DeviceType;
using IndexType = typename Matrix::IndexType;
+ using DiagonalsShiftsType = typename Matrix::DiagonalsShiftsType;
/*
* 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 /
+ * / 1 0 1 0 1 0 \
+ * | 0 1 0 1 0 1 |
+ * | 0 0 1 0 1 0 |
+ * | 1 0 0 1 0 1 |
+ * | 0 1 0 0 1 0 |
+ * | 0 0 1 0 0 1 |
+ * \ 0 0 0 1 0 0 /
*/
const IndexType rows = 7;
const IndexType cols = 6;
- Matrix m( rows, cols );
-
- 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 );
+ Matrix m( rows, cols, DiagonalsShiftsType( { -3, 0, 2, 4 } ) );
+ m.setValue( 1.0 );
+
+ EXPECT_EQ( m.getElement( 0, 0 ), 1 );
+ EXPECT_EQ( m.getElement( 0, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 0, 2 ), 1 );
+ EXPECT_EQ( m.getElement( 0, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 0, 4 ), 1 );
+ EXPECT_EQ( m.getElement( 0, 5 ), 0 );
+
+ EXPECT_EQ( m.getElement( 1, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 1, 1 ), 1 );
+ EXPECT_EQ( m.getElement( 1, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 1, 3 ), 1 );
+ EXPECT_EQ( m.getElement( 1, 4 ), 0 );
+ EXPECT_EQ( m.getElement( 1, 5 ), 1 );
+
+ EXPECT_EQ( m.getElement( 2, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 2, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 2, 2 ), 1 );
+ EXPECT_EQ( m.getElement( 2, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 2, 4 ), 1 );
+ EXPECT_EQ( m.getElement( 2, 5 ), 0 );
+
+ EXPECT_EQ( m.getElement( 3, 0 ), 1 );
+ EXPECT_EQ( m.getElement( 3, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 3, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 3, 3 ), 1 );
+ EXPECT_EQ( m.getElement( 3, 4 ), 0 );
+ EXPECT_EQ( m.getElement( 3, 5 ), 1 );
+
+ EXPECT_EQ( m.getElement( 4, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 1 ), 1 );
+ EXPECT_EQ( m.getElement( 4, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 4 ), 1 );
+ EXPECT_EQ( m.getElement( 4, 5 ), 0 );
+
+ EXPECT_EQ( m.getElement( 5, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 2 ), 1 );
+ EXPECT_EQ( m.getElement( 5, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 4 ), 0 );
+ EXPECT_EQ( m.getElement( 5, 5 ), 1 );
+
+ 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 ), 1 );
+ EXPECT_EQ( m.getElement( 6, 4 ), 0 );
+ EXPECT_EQ( m.getElement( 6, 5 ), 0 );
}
template< typename Matrix >
@@ -429,61 +335,70 @@ void test_SetElement()
using RealType = typename Matrix::RealType;
using DeviceType = typename Matrix::DeviceType;
using IndexType = typename Matrix::IndexType;
+ using DiagonalsShiftsType = typename Matrix::DiagonalsShiftsType;
/*
* 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 /
+ * / 1 2 0 0 5 \
+ * | 0 7 8 0 0 |
+ * | 0 0 13 14 0 |
+ * | 16 0 0 19 20 |
+ * \ 0 22 0 0 25 /
*/
const IndexType rows = 5;
const IndexType cols = 5;
-
- Matrix m( rows, cols );
+ DiagonalsShiftsType diagonals{-3, 0, 1, 4 };
+ Matrix m( rows, cols, diagonals );
RealType value = 1;
for( IndexType i = 0; i < rows; i++ )
for( IndexType j = 0; j < cols; j++ )
{
- if( abs( i - j ) > 1 )
+ bool found( false );
+ for( IndexType k = 0; k < diagonals.getSize(); k++ )
+ {
+ if( i + diagonals[ k ] == j )
+ {
+ m.setElement( i, j, value++ );
+ found = true;
+ break;
+ }
+ }
+ if( ! found )
{
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( 0, 4 ), 5 );
- EXPECT_EQ( m.getElement( 1, 0 ), 6 );
+ EXPECT_EQ( m.getElement( 1, 0 ), 0 );
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, 1 ), 0 );
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, 0 ), 16 );
EXPECT_EQ( m.getElement( 3, 1 ), 0 );
- EXPECT_EQ( m.getElement( 3, 2 ), 18 );
+ EXPECT_EQ( m.getElement( 3, 2 ), 0 );
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, 1 ), 22 );
EXPECT_EQ( m.getElement( 4, 2 ), 0 );
- EXPECT_EQ( m.getElement( 4, 3 ), 24 );
+ EXPECT_EQ( m.getElement( 4, 3 ), 0 );
EXPECT_EQ( m.getElement( 4, 4 ), 25 );
}
@@ -493,123 +408,137 @@ void test_AddElement()
using RealType = typename Matrix::RealType;
using DeviceType = typename Matrix::DeviceType;
using IndexType = typename Matrix::IndexType;
+ using DiagonalsShiftsType = typename Matrix::DiagonalsShiftsType;
/*
- * Sets up the following 6x5 matrix:
+ * 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 |
- * \ 0 0 0 0 30 /
+ * / 1 2 0 0 5 \
+ * | 0 7 8 0 0 |
+ * | 0 0 13 14 0 |
+ * | 0 0 0 19 20 |
+ * \ 0 0 0 0 25 /
*/
-
- const IndexType rows = 6;
+ const IndexType rows = 5;
const IndexType cols = 5;
+ DiagonalsShiftsType diagonals{-3, 0, 1, 4 };
+ Matrix m( rows, cols, diagonals );
- 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 )
- m.setElement( i, j, value );
- value++;
- }
+ RealType value = 1;
+ for( IndexType i = 0; i < rows; i++ )
+ for( IndexType j = 0; j < cols; j++ )
+ {
+ bool found( false );
+ for( IndexType k = 0; k < diagonals.getSize(); k++ )
+ {
+ if( i + diagonals[ k ] == j )
+ {
+ if( j >= i )
+ m.setElement( i, j, value++ );
+ else value++;
+ found = true;
+ break;
+ }
+ }
+ if( ! found )
+ {
+ EXPECT_THROW( m.setElement( i, j, value++ ), std::logic_error );
+ }
+ }
// 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( 0, 4 ), 5 );
- EXPECT_EQ( m.getElement( 1, 0 ), 6 );
+ EXPECT_EQ( m.getElement( 1, 0 ), 0 );
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, 1 ), 0 );
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, 2 ), 0 );
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, 3 ), 0 );
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 /
+ * / 1 2 0 0 5 \ / 1 2 0 0 5 \ / 3 6 0 0 15 \
+ * | 0 7 8 0 0 | | 0 7 8 0 0 | | 0 21 24 0 0 |
+ * 2 * | 0 0 13 14 0 | + | 0 0 13 14 0 | = | 0 0 39 42 0 |
+ * | 0 0 0 19 20 | | 16 0 0 19 20 | | 16 0 0 57 60 |
+ * \ 0 0 0 0 25 / \ 0 22 0 0 25 / \ 0 22 0 0 75 /
+ *
*/
- RealType newValue = 1;
+ value = 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 );
+ {
+ bool found( false );
+ for( IndexType k = 0; k < diagonals.getSize(); k++ )
+ {
+ if( i + diagonals[ k ] == j )
+ {
+ m.addElement( i, j, value++, multiplicator );
+ found = true;
+ break;
+ }
+ }
+ if( ! found )
+ {
+ EXPECT_THROW( m.addElement( i, j, value++, multiplicator ), std::logic_error );
+ }
+ }
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( 0, 4 ), 15 );
- 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, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 1, 1 ), 21 );
+ EXPECT_EQ( m.getElement( 1, 2 ), 24 );
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, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 2, 2 ), 39 );
+ EXPECT_EQ( m.getElement( 2, 3 ), 42 );
EXPECT_EQ( m.getElement( 2, 4 ), 0 );
- EXPECT_EQ( m.getElement( 3, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 3, 0 ), 16 );
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( 3, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 3, 3 ), 57 );
+ EXPECT_EQ( m.getElement( 3, 4 ), 60 );
EXPECT_EQ( m.getElement( 4, 0 ), 0 );
- EXPECT_EQ( m.getElement( 4, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 1 ), 22 );
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 );
+ EXPECT_EQ( m.getElement( 4, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 4 ), 75 );
}
template< typename Matrix >
@@ -618,58 +547,75 @@ void test_SetRow()
using RealType = typename Matrix::RealType;
using DeviceType = typename Matrix::DeviceType;
using IndexType = typename Matrix::IndexType;
+ using DiagonalsShiftsType = typename Matrix::DiagonalsShiftsType;
/*
- * Sets up the following 3x7 matrix:
+ * Sets up the following 5x7 matrix:
*
- * / 1 2 0 0 0 0 0 \
- * | 8 9 10 0 0 0 0 |
- * \ 0 16 17 18 0 0 0 /
+ * / 1 0 2 0 3 0 0 \
+ * | 4 5 0 6 0 7 0 |
+ * | 0 8 9 0 10 0 11 |
+ * | 0 0 12 13 0 14 0 |
+ * \ 0 0 0 15 16 0 17 /
*/
- const IndexType rows = 3;
+ const IndexType rows = 5;
const IndexType cols = 7;
- Matrix m( rows, cols );
+ Matrix m( rows, cols, DiagonalsShiftsType({ -1, 0, 2, 4 }) );
auto matrix_view = m.getView();
auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
- RealType values[ 3 ][ 3 ] {
- { 1, 2, 0 },
- { 8, 9, 10 },
- { 16, 17, 18 } };
+ RealType values[ 5 ][ 4 ] {
+ { 0, 1, 2, 3 },
+ { 4, 5, 6, 7 },
+ { 8, 9, 10, 11 },
+ { 12, 13, 14, 0 },
+ { 15, 16, 17, 0 } };
auto row = matrix_view.getRow( rowIdx );
- for( IndexType i = 0; i < 3; i++ )
- {
- if( rowIdx == 0 && i > 1 )
- break;
+ for( IndexType i = 0; i < 4; i++ )
row.setElement( i, values[ rowIdx ][ i ] );
- }
};
- TNL::Algorithms::ParallelFor< DeviceType >::exec( 0, 3, f );
+ TNL::Algorithms::ParallelFor< DeviceType >::exec( ( IndexType) 0, rows, f );
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, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 0, 2 ), 2 );
EXPECT_EQ( m.getElement( 0, 3 ), 0 );
- EXPECT_EQ( m.getElement( 0, 4 ), 0 );
+ EXPECT_EQ( m.getElement( 0, 4 ), 3 );
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, 0 ), 4 );
+ EXPECT_EQ( m.getElement( 1, 1 ), 5 );
+ EXPECT_EQ( m.getElement( 1, 2 ), 0 );
+ EXPECT_EQ( m.getElement( 1, 3 ), 6 );
EXPECT_EQ( m.getElement( 1, 4 ), 0 );
- EXPECT_EQ( m.getElement( 1, 5 ), 0 );
+ EXPECT_EQ( m.getElement( 1, 5 ), 7 );
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 ), 17 );
- EXPECT_EQ( m.getElement( 2, 3 ), 18 );
- EXPECT_EQ( m.getElement( 2, 4 ), 0 );
+ EXPECT_EQ( m.getElement( 2, 1 ), 8 );
+ EXPECT_EQ( m.getElement( 2, 2 ), 9 );
+ EXPECT_EQ( m.getElement( 2, 3 ), 0 );
+ EXPECT_EQ( m.getElement( 2, 4 ), 10 );
EXPECT_EQ( m.getElement( 2, 5 ), 0 );
- EXPECT_EQ( m.getElement( 2, 6 ), 0 );
+ EXPECT_EQ( m.getElement( 2, 6 ), 11 );
+
+ EXPECT_EQ( m.getElement( 3, 0 ), 0 );
+ EXPECT_EQ( m.getElement( 3, 1 ), 0 );
+ EXPECT_EQ( m.getElement( 3, 2 ), 12 );
+ EXPECT_EQ( m.getElement( 3, 3 ), 13 );
+ EXPECT_EQ( m.getElement( 3, 4 ), 0 );
+ EXPECT_EQ( m.getElement( 3, 5 ), 14 );
+ EXPECT_EQ( m.getElement( 3, 6 ), 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 ), 15 );
+ EXPECT_EQ( m.getElement( 4, 4 ), 16 );
+ EXPECT_EQ( m.getElement( 4, 5 ), 0 );
+ EXPECT_EQ( m.getElement( 4, 6 ), 17 );
}
template< typename Matrix >
@@ -678,27 +624,31 @@ void test_AddRow()
using RealType = typename Matrix::RealType;
using DeviceType = typename Matrix::DeviceType;
using IndexType = typename Matrix::IndexType;
+ using DiagonalsShiftsType = typename Matrix::DiagonalsShiftsType;
+
/*
* 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 /
+ * / 1 2 3 0 0 \
+ * | 0 7 8 9 0 |
+ * | 0 0 13 14 15 |
+ * | 0 0 0 19 20 |
+ * | 0 0 0 0 25 |
+ * \ 0 0 0 0 0 /
*/
const IndexType rows = 6;
const IndexType cols = 5;
+ DiagonalsShiftsType diagonals( { -2, 0, 1, 2 } );
- Matrix m( rows, cols );
+ Matrix m( rows, cols, diagonals );
RealType value = 1;
for( IndexType i = 0; i < rows; i++ )
for( IndexType j = 0; j < cols; j++ )
{
- if( abs( i - j ) <= 1 )
+ IndexType offset = j - i;
+ if( diagonals.containsValue( offset ) && offset >= 0)
m.setElement( i, j, value );
value++;
}
@@ -706,63 +656,63 @@ void test_AddRow()
// 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, 2 ), 3 );
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, 0 ), 0 );
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, 3 ), 9 );
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, 1 ), 0 );
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( 2, 4 ), 15 );
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, 2 ), 0 );
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, 3 ), 0 );
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 );
+ EXPECT_EQ( m.getElement( 5, 4 ), 0 );
// 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:
*
- * / 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 /
+ * / 0 0 0 0 0 0 \ / 1 2 3 0 0 \ / 11 0 0 0 0 \ / 11 0 0 0 0 \
+ * | 0 1 0 0 0 0 | | 0 7 8 9 0 | | 0 22 0 0 0 | | 0 29 8 9 0 |
+ * | 0 0 2 0 0 0 | * | 0 0 13 14 15 | + | 33 0 33 0 0 | = | 33 0 59 28 30 |
+ * | 0 0 0 3 0 0 | | 0 0 0 19 20 | | 0 44 0 44 0 | | 0 44 0 101 60 |
+ * | 0 0 0 0 4 0 | | 0 0 0 0 25 | | 0 0 55 0 55 | | 0 0 55 0 155 |
+ * \ 0 0 0 0 0 5 / \ 0 0 0 0 0 / \ 0 0 0 66 0 / \ 0 0 0 66 0 /
*/
auto matrix_view = m.getView();
auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
- RealType values[ 6 ][ 3 ] {
- { 11, 11, 0 },
- { 22, 22, 22 },
- { 33, 33, 33 },
- { 44, 44, 44 },
- { 55, 55, 55 },
- { 66, 66, 66 } };
+ RealType values[ 6 ][ 4 ] {
+ { 0, 11, 0, 0 },
+ { 0, 22, 0, 0 },
+ { 33, 33, 0, 0 },
+ { 44, 44, 0, 0 },
+ { 55, 55, 0, 0 },
+ { 66, 0, 0, 0 } };
auto row = matrix_view.getRow( rowIdx );
- for( IndexType i = 0; i < 3; i++ )
+ for( IndexType i = 0; i < 4; i++ )
{
RealType& val = row.getValue( i );
val = rowIdx * val + values[ rowIdx ][ i ];
@@ -770,42 +720,41 @@ 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, 1 ), 0 );
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, 0 ), 0 );
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, 2 ), 8 );
+ EXPECT_EQ( m.getElement( 1, 3 ), 9 );
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, 0 ), 33 );
+ EXPECT_EQ( m.getElement( 2, 1 ), 0 );
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( 2, 3 ), 28 );
+ EXPECT_EQ( m.getElement( 2, 4 ), 30 );
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, 1 ), 44 );
+ EXPECT_EQ( m.getElement( 3, 2 ), 0 );
EXPECT_EQ( m.getElement( 3, 3 ), 101 );
- EXPECT_EQ( m.getElement( 3, 4 ), 104 );
+ EXPECT_EQ( m.getElement( 3, 4 ), 60 );
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, 2 ), 55 );
+ EXPECT_EQ( m.getElement( 4, 3 ), 0 );
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 );
+ EXPECT_EQ( m.getElement( 5, 3 ), 66 );
+ EXPECT_EQ( m.getElement( 5, 4 ), 0 );
}
template< typename Matrix >
@@ -814,6 +763,7 @@ void test_VectorProduct()
using RealType = typename Matrix::RealType;
using DeviceType = typename Matrix::DeviceType;
using IndexType = typename Matrix::IndexType;
+ using DiagonalsShiftsType = typename Matrix::DiagonalsShiftsType;
/*
* Sets up the following 5x4 matrix:
@@ -826,8 +776,9 @@ void test_VectorProduct()
*/
const IndexType rows = 5;
const IndexType cols = 4;
+ DiagonalsShiftsType diagonals{ -1, 0, 1 };
- Matrix m( rows, cols );
+ Matrix m( rows, cols, diagonals );
RealType value = 1;
for( IndexType i = 0; i < rows; i++ )
@@ -861,6 +812,8 @@ void test_AddMatrix()
using RealType = typename Matrix1::RealType;
using DeviceType = typename Matrix1::DeviceType;
using IndexType = typename Matrix1::IndexType;
+ using DiagonalsShiftsType1 = typename Matrix1::DiagonalsShiftsType;
+ using DiagonalsShiftsType2 = typename Matrix2::DiagonalsShiftsType;
/*
* Sets up the following 5x4 matrix:
@@ -873,8 +826,10 @@ void test_AddMatrix()
*/
const IndexType rows = 5;
const IndexType cols = 4;
+ DiagonalsShiftsType1 diagonals1;
+ DiagonalsShiftsType2 diagonals2;
- Matrix1 m( rows, cols );
+ Matrix1 m( rows, cols, diagonals1 );
RealType value = 1;
for( IndexType i = 0; i < rows; i++ )
@@ -894,7 +849,7 @@ void test_AddMatrix()
* | 0 0 9 10 |
* \ 0 0 0 11 /
*/
- Matrix2 m2( rows, cols );
+ Matrix2 m2( rows, cols, diagonals2 );
RealType newValue = 1;
for( IndexType i = 0; i < rows; i++ )
@@ -1457,7 +1412,7 @@ TYPED_TEST( MatrixTest, getNonemptyRowsCountTest )
test_GetNonemptyRowsCount< MatrixType >();
}
-/*
+
TYPED_TEST( MatrixTest, getCompressedRowLengthTest )
{
using MatrixType = typename TestFixture::MatrixType;
@@ -1465,13 +1420,6 @@ TYPED_TEST( MatrixTest, getCompressedRowLengthTest )
test_GetCompressedRowLengths< MatrixType >();
}
-TYPED_TEST( MatrixTest, getRowLengthTest )
-{
- using MatrixType = typename TestFixture::MatrixType;
-
- test_GetRowLength< MatrixType >();
-}
-
TYPED_TEST( MatrixTest, getAllocatedElementsCountTest )
{
using MatrixType = typename TestFixture::MatrixType;
@@ -1535,7 +1483,7 @@ TYPED_TEST( MatrixTest, vectorProductTest )
test_VectorProduct< MatrixType >();
}
-TYPED_TEST( MatrixTest, addMatrixTest )
+/*TYPED_TEST( MatrixTest, addMatrixTest )
{
using MatrixType = typename TestFixture::MatrixType;
--
GitLab