Reimplementing tridiagonal matrix. (9a565126) · Commits · TNL / tnl-dev

src/TNL/Matrices/SparseMatrix.hpp

+0 −1

Original line number	Diff line number	Diff line
		@@ -515,7 +515,6 @@ forRows( IndexType first, IndexType last, Function& function )
		return true;
		};
		this->segments.forSegments( first, last, f );

		}

		template< typename Real,

src/TNL/Matrices/Tridiagonal.h

+29 −33

Original line number	Diff line number	Diff line
		@@ -12,15 +12,14 @@

		#include <TNL/Matrices/Matrix.h>
		#include <TNL/Containers/Vector.h>
		#include <TNL/Matrices/TridiagonalRow.h>
		#include <TNL/Matrices/TridiagonalMatrixRowView.h>
		#include <TNL/Containers/Segments/Ellpack.h>
		#include <TNL/Matrices/details/TridiagonalMatrixIndexer.h>
		#include <TNL/Matrices/TridiagonalMatrixView.h>

		namespace TNL {
		namespace Matrices {

		template< typename Device >
		class TridiagonalDeviceDependentCode;

		template< typename Real = double,
		typename Device = Devices::Host,
		typename Index = int,
		@@ -28,27 +27,23 @@ template< typename Real = double,
		typename RealAllocator = typename Allocators::Default< Device >::template Allocator< Real > >
		class Tridiagonal : public Matrix< Real, Device, Index, RealAllocator >
		{
		private:
		// convenient template alias for controlling the selection of copy-assignment operator
		template< typename Device2 >
		using Enabler = std::enable_if< ! std::is_same< Device2, Device >::value >;

		// friend class will be needed for templated assignment operators
		template< typename Real2, typename Device2, typename Index2 >
		friend class Tridiagonal;

		public:
		using RealType = Real;
		using DeviceType = Device;
		using IndexType = Index;
		using RealAllocatorType = RealAllocator;
		using BaseType = Matrix< Real, Device, Index, RealAllocator >;
		using IndexerType = details::TridiagonalMatrixIndexer< IndexType, RowMajorOrder >;
		using ValuesType = typename BaseType::ValuesVector;
		using ValuesViewType = typename ValuesType::ViewType;
		//using ViewType = TridiagonalMatrixView< Real, Device, Index, RowMajorOrder >;
		//using ConstViewType = TridiagonalMatrixView< typename std::add_const< Real >::type, Device, Index, RowMajorOrder >;
		using RowView = TridiagonalMatrixRowView< SegmentViewType, ValuesViewType >;
		using ViewType = TridiagonalMatrixView< Real, Device, Index, RowMajorOrder >;
		using ConstViewType = TridiagonalMatrixView< typename std::add_const< Real >::type, Device, Index, RowMajorOrder >;
		using RowView = TridiagonalMatrixRowView< 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,
		@@ -70,7 +65,8 @@ class Tridiagonal : public Matrix< Real, Device, Index, RealAllocator >
		void setDimensions( const IndexType rows,
		const IndexType columns );

		void setCompressedRowLengths( ConstCompressedRowLengthsVectorView rowLengths );
		//template< typename Vector >
		void setCompressedRowLengths( const ConstCompressedRowLengthsVectorView rowCapacities );

		template< typename Vector >
		void getCompressedRowLengths( Vector& rowLengths ) const;
		@@ -80,8 +76,8 @@ class Tridiagonal : public Matrix< Real, Device, Index, RealAllocator >

		IndexType getMaxRowLength() const;

		template< typename Real2, typename Device2, typename Index2 >
		void setLike( const Tridiagonal< Real2, Device2, Index2 >& m );
		template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_, typename RealAllocator_ >
		void setLike( const Tridiagonal< Real_, Device_, Index_, RowMajorOrder_, RealAllocator_ >& m );

		IndexType getNumberOfMatrixElements() const;

		@@ -91,11 +87,15 @@ class Tridiagonal : public Matrix< Real, Device, Index, RealAllocator >

		void reset();

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

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

		template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_ >
		bool operator != ( const Tridiagonal< Real_, Device_, Index_ >& matrix ) const;
		RowView getRow( const IndexType& rowIdx );

		const RowView getRow( const IndexType& rowIdx ) const;

		void setValue( const RealType& v );

		@@ -139,8 +139,8 @@ class Tridiagonal : public Matrix< Real, Device, Index, RealAllocator >
		void vectorProduct( const InVector& inVector,
		OutVector& outVector ) const;

		template< typename Real2, typename Index2 >
		void addMatrix( const Tridiagonal< Real2, Device, Index2 >& matrix,
		template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_, typename RealAllocator_ >
		void addMatrix( const Tridiagonal< Real_, Device_, Index_, RowMajorOrder_, RealAllocator_ >& matrix,
		const RealType& matrixMultiplicator = 1.0,
		const RealType& thisMatrixMultiplicator = 1.0 );

		@@ -159,9 +159,8 @@ class Tridiagonal : public Matrix< Real, Device, Index, RealAllocator >
		Tridiagonal& operator=( const Tridiagonal& matrix );

		// cross-device copy assignment
		template< typename Real2, typename Device2, typename Index2,
		typename = typename Enabler< Device2 >::type >
		Tridiagonal& operator=( const Tridiagonal< Real2, Device2, Index2 >& matrix );
		template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_, typename RealAllocator_ >
		Tridiagonal& operator=( const Tridiagonal< Real_, Device_, Index_, RowMajorOrder_, RealAllocator_ >& matrix );

		void save( File& file ) const;

		@@ -177,12 +176,9 @@ class Tridiagonal : public Matrix< Real, Device, Index, RealAllocator >

		__cuda_callable__
		IndexType getElementIndex( const IndexType row,
		const IndexType column ) const;

		Containers::Vector< RealType, DeviceType, IndexType > values;
		const IndexType localIdx ) const;

		typedef TridiagonalDeviceDependentCode< DeviceType > DeviceDependentCode;
		friend class TridiagonalDeviceDependentCode< DeviceType >;
		IndexerType indexer;
		};

		} // namespace Matrices

src/TNL/Matrices/Tridiagonal.hpp

+228 −75

Original line number	Diff line number	Diff line
		/***************************************************************************
		Tridiagonal_impl.h - description
		Tridiagonal.hpp - description
		-------------------
		begin : Nov 30, 2013
		copyright : (C) 2013 by Tomas Oberhuber
		@@ -41,6 +41,30 @@ Tridiagonal( const IndexType rows, const IndexType columns )
		this->setDimensions( rows, columns );
		}

		template< typename Real,
		typename Device,
		typename Index,
		bool RowMajorOrder,
		typename RealAllocator >
		auto
		Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
		getView() -> ViewType
		{
		return ViewType( this->values.getView(), indexer );
		}

		template< typename Real,
		typename Device,
		typename Index,
		bool RowMajorOrder,
		typename RealAllocator >
		auto
		Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
		getConstView() const -> ConstViewType
		{
		return ConstViewType( this->values.getConstView(), indexer );
		}

		template< typename Real,
		typename Device,
		typename Index,
		@@ -78,7 +102,8 @@ Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
		setDimensions( const IndexType rows, const IndexType columns )
		{
		Matrix< Real, Device, Index >::setDimensions( rows, columns );
		values.setSize( 3*min( rows, columns ) );
		this->indexer.setDimensions( rows, columns );
		this->values.setSize( this->indexer.getStorageSize() );
		this->values = 0.0;
		}

		@@ -87,24 +112,24 @@ template< typename Real,
		typename Index,
		bool RowMajorOrder,
		typename RealAllocator >
		// template< typename Vector >
		void
		Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
		setCompressedRowLengths( ConstCompressedRowLengthsVectorView rowLengths )
		setCompressedRowLengths( const ConstCompressedRowLengthsVectorView rowLengths )
		{
		if( rowLengths[ 0 ] > 2 )
		if( max( rowLengths ) > 3 )
		throw std::logic_error( "Too many non-zero elements per row in a tri-diagonal matrix." );
		const IndexType diagonalLength = min( this->getRows(), this->getColumns() );
		for( Index i = 1; i < diagonalLength-1; i++ )
		if( rowLengths[ i ] > 3 )
		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[ this->getRows()-1 ] > 1 )
		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[ this->getRows()-1 ] > 2 )
		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[ this->getRows()-1 ] > 3 )
		if( rowLengths.getElement( this->getRows()-1 ) > 3 )
		throw std::logic_error( "Too many non-zero elements per row in a tri-diagonal matrix." );
		}

		@@ -146,10 +171,10 @@ template< typename Real,
		typename Index,
		bool RowMajorOrder,
		typename RealAllocator >
		template< typename Real2, typename Device2, typename Index2 >
		template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_, typename RealAllocator_ >
		void
		Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
		setLike( const Tridiagonal< Real2, Device2, Index2 >& m )
		setLike( const Tridiagonal< Real_, Device_, Index_, RowMajorOrder_, RealAllocator_ >& m )
		{
		this->setDimensions( m.getRows(), m.getColumns() );
		}
		@@ -250,6 +275,32 @@ setValue( const RealType& v )
		this->values = v;
		}

		template< typename Real,
		typename Device,
		typename Index,
		bool RowMajorOrder,
		typename RealAllocator >
		__cuda_callable__
		auto
		Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
		getRow( const IndexType& rowIdx ) const -> const RowView
		{
		return RowView( this->values.getView(), this->indexer );
		}

		template< typename Real,
		typename Device,
		typename Index,
		bool RowMajorOrder,
		typename RealAllocator >
		__cuda_callable__
		auto
		Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
		getRow( const IndexType& rowIdx ) -> RowView
		{
		return RowView( this->values.getView(), this->indexer );
		}

		template< typename Real,
		typename Device,
		typename Index,
		@@ -259,6 +310,12 @@ bool
		Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
		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 )
		throw std::logic_error( "Wrong matrix element coordinates in tridiagonal matrix." );
		this->values.setElement( this->getElementIndex( row, column ), value );
		return true;
		}
		@@ -275,6 +332,12 @@ addElement( const IndexType row,
		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 )
		throw std::logic_error( "Wrong matrix element coordinates in tridiagonal matrix." );
		const Index i = this->getElementIndex( row, column );
		this->values.setElement( i, thisElementMultiplicator * this->values.getElement( i ) + value );
		return true;
		@@ -289,6 +352,11 @@ Real
		Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
		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 ) );
		@@ -304,7 +372,46 @@ void
		Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
		rowsReduction( IndexType first, IndexType last, Fetch& fetch, Reduce& reduce, Keep& keep, const FetchReal& zero ) const
		{

		const auto values_view = this->values.getConstView();
		const auto indexer_ = this->indexer;
		const auto rows = this->getRows();
		const auto columns = this->getColumns();
		const auto size = this->size;
		auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
		//bool compute;
		if( rowIdx == 0 )
		{
		IndexType i_0 = indexer.getGlobalIndex( 0, 0 );
		IndexType i_1 = indexer.getGlobalIndex( 0, 1 );
		keep( 0, reduce( fetch( 0, 0, i_0, values_view[ i_0 ] ),
		fetch( 0, 1, i_1, values_view[ i_1 ] ) ) );
		return;
		}
		if( rowIdx < size \|\| columns > rows )
		{
		IndexType i_0 = indexer.getGlobalIndex( rowIdx, 0 );
		IndexType i_1 = indexer.getGlobalIndex( rowIdx, 1 );
		IndexType i_2 = indexer.getGlobalIndex( rowIdx, 2 );

		keep( rowIdx, reduce( reduce( fetch( rowIdx, rowIdx - 1, i_0, values_view[ i_0 ] ),
		fetch( rowIdx, rowIdx, i_1, values_view[ i_1 ] ) ),
		fetch( rowIdx, rowIdx + 1, i_2, values_view[ i_2] ) ) );
		return;
		}
		if( rows == columns )
		{
		IndexType i_0 = indexer.getGlobalIndex( rowIdx, 0 );
		IndexType i_1 = indexer.getGlobalIndex( rowIdx, 1 );
		keep( rowIdx, reduce( fetch( rowIdx, rowIdx - 1, i_0, values_view[ i_0 ] ),
		fetch( rowIdx, rowIdx, i_1, values_view[ i_1 ] ) ) );
		}
		else
		{
		IndexType i_0 = indexer.getGlobalIndex( rowIdx, 0 );
		keep( rowIdx, fetch( rowIdx, rowIdx, i_0, values_view[ i_0 ] ) );
		}
		};
		Algorithms::ParallelFor< DeviceType >::exec( first, last, f );
		}

		template< typename Real,
		@@ -317,7 +424,7 @@ void
		Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
		allRowsReduction( Fetch& fetch, Reduce& reduce, Keep& keep, const FetchReal& zero ) const
		{

		this->rowsReduction( 0, this->getRows(), fetch, reduce, keep, zero );
		}

		template< typename Real,
		@@ -330,7 +437,45 @@ void
		Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
		forRows( IndexType first, IndexType last, Function& function ) const
		{

		const auto values_view = this->values.getConstView();
		const auto indexer_ = this->indexer;
		const auto rows = this->getRows();
		const auto columns = this->getColumns();
		const auto size = this->size;
		auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
		//bool compute;
		if( rowIdx == 0 )
		{
		IndexType i_0 = indexer.getGlobalIndex( 0, 0 );
		IndexType i_1 = indexer.getGlobalIndex( 0, 1 );
		function( 0, 1, rowIdx, values_view[ i_0 ] );
		function( 0, 2, rowIdx + 1, values_view[ i_1 ] );
		return;
		}
		if( rowIdx < size \|\| columns > rows )
		{
		IndexType i_0 = indexer.getGlobalIndex( rowIdx, 0 );
		IndexType i_1 = indexer.getGlobalIndex( rowIdx, 1 );
		IndexType i_2 = indexer.getGlobalIndex( rowIdx, 2 );
		function( rowIdx, 0, rowIdx - 1, values_view[ i_0 ] );
		function( rowIdx, 1, rowIdx, values_view[ i_1 ] );
		function( rowIdx, 2, rowIdx + 1, values_view[ i_2 ] );
		return;
		}
		if( rows == columns )
		{
		IndexType i_0 = indexer.getGlobalIndex( rowIdx, 0 );
		IndexType i_1 = indexer.getGlobalIndex( rowIdx, 1 );
		function( rowIdx, 0, rowIdx - 1, values_view[ i_0 ] );
		function( rowIdx, 1, rowIdx, values_view[ i_1 ] );
		}
		else
		{
		IndexType i_0 = indexer.getGlobalIndex( rowIdx, 0 );
		function( rowIdx, 0, rowIdx, values_view[ i_0 ] );
		}
		};
		Algorithms::ParallelFor< DeviceType >::exec( first, last, f );
		}

		template< typename Real,
		@@ -343,6 +488,45 @@ void
		Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
		forRows( IndexType first, IndexType last, Function& function )
		{
		const auto values_view = this->values.getConstView();
		const auto indexer_ = this->indexer;
		const auto rows = this->getRows();
		const auto columns = this->getColumns();
		const auto size = this->size;
		auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
		//bool compute;
		if( rowIdx == 0 )
		{
		IndexType i_0 = indexer.getGlobalIndex( 0, 0 );
		IndexType i_1 = indexer.getGlobalIndex( 0, 1 );
		function( 0, 1, rowIdx, values_view[ i_0 ] );
		function( 0, 2, rowIdx + 1, values_view[ i_1 ] );
		return;
		}
		if( rowIdx < size \|\| columns > rows )
		{
		IndexType i_0 = indexer.getGlobalIndex( rowIdx, 0 );
		IndexType i_1 = indexer.getGlobalIndex( rowIdx, 1 );
		IndexType i_2 = indexer.getGlobalIndex( rowIdx, 2 );
		function( rowIdx, 0, rowIdx - 1, values_view[ i_0 ] );
		function( rowIdx, 1, rowIdx, values_view[ i_1 ] );
		function( rowIdx, 2, rowIdx + 1, values_view[ i_2 ] );
		return;
		}
		if( rows == columns )
		{
		IndexType i_0 = indexer.getGlobalIndex( rowIdx, 0 );
		IndexType i_1 = indexer.getGlobalIndex( rowIdx, 1 );
		function( rowIdx, 0, rowIdx - 1, values_view[ i_0 ] );
		function( rowIdx, 1, rowIdx, values_view[ i_1 ] );
		}
		else
		{
		IndexType i_0 = indexer.getGlobalIndex( rowIdx, 0 );
		function( rowIdx, 0, rowIdx, values_view[ i_0 ] );
		}
		};
		Algorithms::ParallelFor< DeviceType >::exec( first, last, f );
		}

		template< typename Real,
		@@ -355,7 +539,7 @@ void
		Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
		forAllRows( Function& function ) const
		{

		this->forRows( 0, this->getRows(), function );
		}

		template< typename Real,
		@@ -368,45 +552,9 @@ void
		Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
		forAllRows( Function& function )
		{

		}

		template< typename Real,
		typename Device,
		typename Index,
		bool RowMajorOrder,
		typename RealAllocator >
		__cuda_callable__
		typename Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::MatrixRow
		Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
		getRow( const IndexType rowIndex )
		{
		if( std::is_same< Device, Devices::Host >::value )
		return MatrixRow( &this->values.getData()[ this->getElementIndex( rowIndex, rowIndex ) ],
		rowIndex,
		this->getColumns(),
		1 );
		if( std::is_same< Device, Devices::Cuda >::value )
		return MatrixRow( &this->values.getData()[ this->getElementIndex( rowIndex, rowIndex ) ],
		rowIndex,
		this->getColumns(),
		this->rows );
		}

		template< typename Real,
		typename Device,
		typename Index,
		bool RowMajorOrder,
		typename RealAllocator >
		__cuda_callable__
		const typename Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::MatrixRow
		Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
		getRow( const IndexType rowIndex ) const
		{
		throw Exceptions::NotImplementedError();
		this->forRows( 0, this->getRows(), function );
		}


		template< typename Real,
		typename Device,
		typename Index,
		@@ -414,8 +562,9 @@ template< typename Real,
		typename RealAllocator >
		template< typename Vector >
		__cuda_callable__
		typename Vector::RealType Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::rowVectorProduct( const IndexType row,
		const Vector& vector ) const
		typename Vector::RealType
		Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
		rowVectorProduct( const IndexType row, const Vector& vector ) const
		{
		return TridiagonalDeviceDependentCode< Device >::
		rowVectorProduct( this->rows,
		@@ -431,8 +580,9 @@ template< typename Real,
		typename RealAllocator >
		template< typename InVector,
		typename OutVector >
		void Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::vectorProduct( const InVector& inVector,
		OutVector& outVector ) const
		void
		Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
		vectorProduct( const InVector& inVector, OutVector& outVector ) const
		{
		TNL_ASSERT( this->getColumns() == inVector.getSize(),
		std::cerr << "Matrix columns: " << this->getColumns() << std::endl
		@@ -441,7 +591,7 @@ void Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::vectorPro
		std::cerr << "Matrix rows: " << this->getRows() << std::endl
		<< "Vector size: " << outVector.getSize() << std::endl );

		DeviceDependentCode::vectorProduct( *this, inVector, outVector );
		//DeviceDependentCode::vectorProduct( *this, inVector, outVector );
		}

		template< typename Real,
		@@ -449,8 +599,10 @@ template< typename Real,
		typename Index,
		bool RowMajorOrder,
		typename RealAllocator >
		template< typename Real2, typename Index2 >
		void Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::addMatrix( const Tridiagonal< Real2, Device, Index2 >& matrix,
		template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_, typename RealAllocator_ >
		void
		Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
		addMatrix( const Tridiagonal< Real_, Device_, Index_, RowMajorOrder_, RealAllocator_ >& matrix,
		const RealType& matrixMultiplicator,
		const RealType& thisMatrixMultiplicator )
		{
		@@ -582,13 +734,14 @@ template< typename Real,
		typename Index,
		bool RowMajorOrder,
		typename RealAllocator >
		template< typename Real2, typename Device2, typename Index2, typename >
		template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_, typename RealAllocator_ >
		Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >&
		Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::operator=( const Tridiagonal< Real2, Device2, Index2 >& matrix )
		Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
		operator=( const Tridiagonal< 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< Device2, Devices::Host >::value \|\| std::is_same< Device2, Devices::Cuda >::value,
		static_assert( std::is_same< Device_, Devices::Host >::value \|\| std::is_same< Device_, Devices::Cuda >::value,
		"unknown device" );

		this->setLike( matrix );
		@@ -605,7 +758,6 @@ template< typename Real,
		void Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::save( File& file ) const
		{
		Matrix< Real, Device, Index >::save( file );
		file << this->values;
		}

		template< typename Real,
		@@ -616,7 +768,7 @@ template< typename Real,
		void Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::load( File& file )
		{
		Matrix< Real, Device, Index >::load( file );
		file >> this->values;
		this->indexer.setDimensions( this->getRows(), this->getColumns() );
		}

		template< typename Real,
		@@ -662,17 +814,17 @@ template< typename Real,
		bool RowMajorOrder,
		typename RealAllocator >
		__cuda_callable__
		Index Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::getElementIndex( const IndexType row,
		const IndexType column ) const
		Index Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
		getElementIndex( const IndexType row, const IndexType localIdx ) const
		{
		TNL_ASSERT( row >= 0 && column >= 0 && row < this->rows && column < this->rows,
		std::cerr << " this->rows = " << this->rows
		<< " row = " << row << " column = " << column );
		TNL_ASSERT( abs( row - column ) < 2,
		std::cerr << "row = " << row << " column = " << column << std::endl );
		return TridiagonalDeviceDependentCode< Device >::getElementIndex( this->rows, row, column );
		TNL_ASSERT_GE( row, 0, "" );
		TNL_ASSERT_LT( row, this->getRows(), "" );
		TNL_ASSERT_GE( localIdx, 0, "" );
		TNL_ASSERT_LT( localIdx, 3, "" );
		return this->indexer.getGlobalIndex( row, localIdx );
		}

		/*
		template<>
		class TridiagonalDeviceDependentCode< Devices::Host >
		{
		@@ -774,6 +926,7 @@ class TridiagonalDeviceDependentCode< Devices::Cuda >
		MatrixVectorProductCuda( matrix, inVector, outVector );
		}
		};
		*/

		} // namespace Matrices
		} // namespace TNL

src/TNL/Matrices/TridiagonalRow.h→src/TNL/Matrices/TridiagonalMatrixRowView.h

+59 −0

File changed and moved.

Preview size limit exceeded, changes collapsed.

src/TNL/Matrices/TridiagonalMatrixRowView.hpp

0 → 100644

+75 −0

File added.

Preview size limit exceeded, changes collapsed.