From f4d7c074e0383c9a9aa908e2de833950ceacd173 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= <oberhuber.tomas@gmail.com>
Date: Thu, 4 Feb 2021 15:15:57 +0100
Subject: [PATCH] Added getRowCapacities method to matrices.

---
 src/TNL/Matrices/DenseMatrix.h               |  10 ++
 src/TNL/Matrices/DenseMatrix.hpp             |  13 ++
 src/TNL/Matrices/DenseMatrixView.h           |  10 ++
 src/TNL/Matrices/DenseMatrixView.hpp         |  14 ++
 src/TNL/Matrices/LambdaMatrix.h              |  10 ++
 src/TNL/Matrices/LambdaMatrix.hpp            |  14 +-
 src/TNL/Matrices/MultidiagonalMatrix.h       |  12 +-
 src/TNL/Matrices/MultidiagonalMatrix.hpp     |  16 +-
 src/TNL/Matrices/MultidiagonalMatrixView.h   |  96 ++++++-----
 src/TNL/Matrices/MultidiagonalMatrixView.hpp |  20 ++-
 src/TNL/Matrices/SparseMatrix.h              |  11 ++
 src/TNL/Matrices/SparseMatrix.hpp            |  32 +++-
 src/TNL/Matrices/SparseMatrixView.h          |  83 +++++-----
 src/TNL/Matrices/SparseMatrixView.hpp        |  23 +++
 src/TNL/Matrices/TridiagonalMatrix.h         |  10 ++
 src/TNL/Matrices/TridiagonalMatrix.hpp       |  13 ++
 src/TNL/Matrices/TridiagonalMatrixView.h     | 163 ++++++++++---------
 src/TNL/Matrices/TridiagonalMatrixView.hpp   |  13 ++
 18 files changed, 395 insertions(+), 168 deletions(-)

diff --git a/src/TNL/Matrices/DenseMatrix.h b/src/TNL/Matrices/DenseMatrix.h
index c12a4347fd..76c0f86252 100644
--- a/src/TNL/Matrices/DenseMatrix.h
+++ b/src/TNL/Matrices/DenseMatrix.h
@@ -225,6 +225,16 @@ class DenseMatrix : public Matrix< Real, Device, Index, RealAllocator >
       template< typename RowCapacitiesVector >
       void setRowCapacities( const RowCapacitiesVector& rowCapacities );
 
+      /**
+       * \brief Compute capacities of all rows.
+       *
+       * The row capacities are not stored explicitly and must be computed.
+       *
+       * \param rowCapacities is a vector where the row capacities will be stored.
+       */
+      template< typename Vector >
+      void getRowCapacities( Vector& rowCapacities ) const;
+
       /**
        * \brief This method recreates the dense matrix from 2D initializer list.
        *
diff --git a/src/TNL/Matrices/DenseMatrix.hpp b/src/TNL/Matrices/DenseMatrix.hpp
index d6d6cb04f3..f0bd56141c 100644
--- a/src/TNL/Matrices/DenseMatrix.hpp
+++ b/src/TNL/Matrices/DenseMatrix.hpp
@@ -192,6 +192,19 @@ setRowCapacities( const RowCapacitiesVector& rowCapacities )
    TNL_ASSERT_LE( max( rowCapacities ), this->getColumns(), "" );
 }
 
+template< typename Real,
+          typename Device,
+          typename Index,
+          ElementsOrganization Organization,
+          typename RealAllocator >
+   template< typename Vector >
+void
+DenseMatrix< Real, Device, Index, Organization, RealAllocator >::
+getRowCapacities( Vector& rowLengths ) const
+{
+   this->view.getCompressedRowLengths( rowLengths );
+}
+
 template< typename Real,
           typename Device,
           typename Index,
diff --git a/src/TNL/Matrices/DenseMatrixView.h b/src/TNL/Matrices/DenseMatrixView.h
index e21d79bdfd..1d54e04f35 100644
--- a/src/TNL/Matrices/DenseMatrixView.h
+++ b/src/TNL/Matrices/DenseMatrixView.h
@@ -174,6 +174,16 @@ class DenseMatrixView : public MatrixView< Real, Device, Index >
        */
       virtual String getSerializationTypeVirtual() const;
 
+      /**
+       * \brief Compute capacities of all rows.
+       *
+       * The row capacities are not stored explicitly and must be computed.
+       *
+       * \param rowCapacities is a vector where the row capacities will be stored.
+       */
+      template< typename Vector >
+      void getRowCapacities( Vector& rowCapacities ) const;
+
       /**
        * \brief Computes number of non-zeros in each row.
        *
diff --git a/src/TNL/Matrices/DenseMatrixView.hpp b/src/TNL/Matrices/DenseMatrixView.hpp
index a19edf6452..36ade91c87 100644
--- a/src/TNL/Matrices/DenseMatrixView.hpp
+++ b/src/TNL/Matrices/DenseMatrixView.hpp
@@ -98,6 +98,20 @@ getSerializationTypeVirtual() const
    return this->getSerializationType();
 }
 
+template< typename Real,
+          typename Device,
+          typename Index,
+          ElementsOrganization Organization >
+   template< typename Vector >
+void
+DenseMatrixView< Real, Device, Index, Organization >::
+getRowCapacities( Vector& rowCapacities ) const
+{
+   rowCapacities.setSize( this->getRows() );
+   rowCapacities = this->getColumns();
+}
+
+
 template< typename Real,
           typename Device,
           typename Index,
diff --git a/src/TNL/Matrices/LambdaMatrix.h b/src/TNL/Matrices/LambdaMatrix.h
index 46f6184b5c..e2a7bc5eb1 100644
--- a/src/TNL/Matrices/LambdaMatrix.h
+++ b/src/TNL/Matrices/LambdaMatrix.h
@@ -148,6 +148,16 @@ class LambdaMatrix
       __cuda_callable__
       IndexType getColumns() const;
 
+      /**
+       * \brief Compute capacities of all rows.
+       *
+       * The row capacities are not stored explicitly and must be computed.
+       *
+       * \param rowCapacities is a vector where the row capacities will be stored.
+       */
+      template< typename Vector >
+      void getRowCapacities( Vector& rowCapacities ) const;
+
       /**
        * \brief Computes number of non-zeros in each row.
        *
diff --git a/src/TNL/Matrices/LambdaMatrix.hpp b/src/TNL/Matrices/LambdaMatrix.hpp
index 65fe9ce93b..3f58446bdd 100644
--- a/src/TNL/Matrices/LambdaMatrix.hpp
+++ b/src/TNL/Matrices/LambdaMatrix.hpp
@@ -87,6 +87,19 @@ getColumns() const
    return this->columns;
 }
 
+template< typename MatrixElementsLambda,
+          typename CompressedRowLengthsLambda,
+          typename Real,
+          typename Device,
+          typename Index >
+   template< typename Vector >
+void
+LambdaMatrix< MatrixElementsLambda, CompressedRowLengthsLambda, Real, Device, Index >::
+getRowCapacities( Vector& rowCapacities ) const
+{
+   this->getCompressedRowLengths( rowCapacities );
+}
+
 template< typename MatrixElementsLambda,
           typename CompressedRowLengthsLambda,
           typename Real,
@@ -344,7 +357,6 @@ performSORIteration( const Vector1& b,
                           Vector2& x,
                           const RealType& omega ) const
 {
-   
 }
 
 template< typename MatrixElementsLambda,
diff --git a/src/TNL/Matrices/MultidiagonalMatrix.h b/src/TNL/Matrices/MultidiagonalMatrix.h
index 341c5c3766..d8f076d000 100644
--- a/src/TNL/Matrices/MultidiagonalMatrix.h
+++ b/src/TNL/Matrices/MultidiagonalMatrix.h
@@ -346,7 +346,7 @@ class MultidiagonalMatrix : public Matrix< Real, Device, Index, RealAllocator >
        *
        * \return Number of diagonals.
        */
-      const IndexType& getDiagonalsCount() const;
+      const IndexType getDiagonalsCount() const;
 
       /**
        * \brief Returns vector with diagonals offsets.
@@ -373,6 +373,16 @@ class MultidiagonalMatrix : public Matrix< Real, Device, Index, RealAllocator >
       template< typename ListReal >
       void setElements( const std::initializer_list< std::initializer_list< ListReal > >& data );
 
+      /**
+       * \brief Compute capacities of all rows.
+       *
+       * The row capacities are not stored explicitly and must be computed.
+       *
+       * \param rowCapacities is a vector where the row capacities will be stored.
+       */
+      template< typename Vector >
+      void getRowCapacities( Vector& rowCapacities ) const;
+
       /**
        * \brief Computes number of non-zeros in each row.
        *
diff --git a/src/TNL/Matrices/MultidiagonalMatrix.hpp b/src/TNL/Matrices/MultidiagonalMatrix.hpp
index 61986d2634..0db276d370 100644
--- a/src/TNL/Matrices/MultidiagonalMatrix.hpp
+++ b/src/TNL/Matrices/MultidiagonalMatrix.hpp
@@ -207,6 +207,20 @@ setRowCapacities( const RowCapacitiesVector& rowLengths )
          throw std::logic_error( "Too many non-zero elements per row in a tri-diagonal matrix." );
 }
 
+template< typename Real,
+          typename Device,
+          typename Index,
+          ElementsOrganization Organization,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Vector >
+void
+MultidiagonalMatrix< Real, Device, Index, Organization, RealAllocator, IndexAllocator >::
+getRowCapacities( Vector& rowCapacities ) const
+{
+   return this->view.getRowCapacities( rowCapacities );
+}
+
 template< typename Real,
           typename Device,
           typename Index,
@@ -248,7 +262,7 @@ template< typename Real,
           ElementsOrganization Organization,
           typename RealAllocator,
           typename IndexAllocator >
-const Index&
+const Index
 MultidiagonalMatrix< Real, Device, Index, Organization, RealAllocator, IndexAllocator >::
 getDiagonalsCount() const
 {
diff --git a/src/TNL/Matrices/MultidiagonalMatrixView.h b/src/TNL/Matrices/MultidiagonalMatrixView.h
index 43882f8266..d2b54ded76 100644
--- a/src/TNL/Matrices/MultidiagonalMatrixView.h
+++ b/src/TNL/Matrices/MultidiagonalMatrixView.h
@@ -163,7 +163,17 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
        * \return Number of diagonals.
        */
       __cuda_callable__
-      const IndexType& getDiagonalsCount() const;
+      const IndexType getDiagonalsCount() const;
+
+      /**
+       * \brief Compute capacities of all rows.
+       *
+       * The row capacities are not stored explicitly and must be computed.
+       *
+       * \param rowCapacities is a vector where the row capacities will be stored.
+       */
+      template< typename Vector >
+      void getRowCapacities( Vector& rowCapacities ) const;
 
       /**
        * \brief Computes number of non-zeros in each row.
@@ -262,14 +272,14 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Set all matrix elements to given value.
-       * 
+       *
        * \param value is the new value of all matrix elements.
        */
       void setValue( const RealType& v );
 
       /**
        * \brief Sets element at given \e row and \e column to given \e value.
-       * 
+       *
        * This method can be called from the host system (CPU) no matter
        * where the matrix is allocated. If the matrix is allocated on GPU this method
        * can be called even from device kernels. If the matrix is allocated in GPU device
@@ -277,11 +287,11 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
        * performance is very low. For higher performance see. \ref MultidiagonalMatrix::getRow
        * or \ref MultidiagonalMatrix::forRows and \ref MultidiagonalMatrix::forAllRows.
        * The call may fail if the matrix row capacity is exhausted.
-       * 
+       *
        * \param row is row index of the element.
        * \param column is columns index of the element.
        * \param value is the value the element will be set to.
-       * 
+       *
        * \par Example
        * \include Matrices/MultidiagonalMatrix/MultidiagonalMatrixViewExample_setElement.cpp
        * \par Output
@@ -294,7 +304,7 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Add element at given \e row and \e column to given \e value.
-       * 
+       *
        * This method can be called from the host system (CPU) no matter
        * where the matrix is allocated. If the matrix is allocated on GPU this method
        * can be called even from device kernels. If the matrix is allocated in GPU device
@@ -302,18 +312,17 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
        * performance is very low. For higher performance see. \ref MultidiagonalMatrix::getRow
        * or \ref MultidiagonalMatrix::forRows and \ref MultidiagonalMatrix::forAllRows.
        * The call may fail if the matrix row capacity is exhausted.
-       * 
+       *
        * \param row is row index of the element.
        * \param column is columns index of the element.
        * \param value is the value the element will be set to.
        * \param thisElementMultiplicator is multiplicator the original matrix element
        *   value is multiplied by before addition of given \e value.
-       * 
+       *
        * \par Example
        * \include Matrices/MultidiagonalMatrix/MultidiagonalMatrixViewExample_addElement.cpp
        * \par Output
        * \include MultidiagonalMatrixViewExample_addElement.out
-       * 
        */
       __cuda_callable__
       void addElement( const IndexType row,
@@ -323,24 +332,23 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Returns value of matrix element at position given by its row and column index.
-       * 
+       *
        * This method can be called from the host system (CPU) no matter
        * where the matrix is allocated. If the matrix is allocated on GPU this method
        * can be called even from device kernels. If the matrix is allocated in GPU device
        * this method is called from CPU, it transfers values of each matrix element separately and so the
        * performance is very low. For higher performance see. \ref MultidiagonalMatrix::getRow
        * or \ref MultidiagonalMatrix::forRows and \ref MultidiagonalMatrix::forAllRows.
-       * 
+       *
        * \param row is a row index of the matrix element.
        * \param column i a column index of the matrix element.
-       * 
+       *
        * \return value of given matrix element.
-       * 
+       *
        * \par Example
        * \include Matrices/MultidiagonalMatrix/MultidiagonalMatrixViewExample_getElement.cpp
        * \par Output
        * \include MultidiagonalMatrixViewExample_getElement.out
-       * 
        */
       __cuda_callable__
       RealType getElement( const IndexType row,
@@ -348,7 +356,7 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Method for performing general reduction on matrix rows for constant instances.
-       * 
+       *
        * \tparam Fetch is a type of lambda function for data fetch declared as
        *          `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`.
        *          The return type of this lambda can be any non void.
@@ -357,14 +365,14 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
        * \tparam Keep is a type of lambda function for storing results of reduction in each row.
        *          It is declared as `keep( const IndexType rowIdx, const double& value )`.
        * \tparam FetchValue is type returned by the Fetch lambda function.
-       * 
+       *
        * \param begin defines beginning of the range [begin,end) of rows to be processed.
        * \param end defines ending of the range [begin,end) of rows to be processed.
        * \param fetch is an instance of lambda function for data fetch.
        * \param reduce is an instance of lambda function for reduction.
        * \param keep in an instance of lambda function for storing results.
        * \param zero is zero of given reduction operation also known as idempotent element.
-       * 
+       *
        * \par Example
        * \include Matrices/MultidiagonalMatrix/MultidiagonalMatrixViewExample_rowsReduction.cpp
        * \par Output
@@ -375,7 +383,7 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Method for performing general reduction on matrix rows.
-       * 
+       *
        * \tparam Fetch is a type of lambda function for data fetch declared as
        *          `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`.
        *          The return type of this lambda can be any non void.
@@ -384,14 +392,14 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
        * \tparam Keep is a type of lambda function for storing results of reduction in each row.
        *          It is declared as `keep( const IndexType rowIdx, const double& value )`.
        * \tparam FetchValue is type returned by the Fetch lambda function.
-       * 
+       *
        * \param begin defines beginning of the range [begin,end) of rows to be processed.
        * \param end defines ending of the range [begin,end) of rows to be processed.
        * \param fetch is an instance of lambda function for data fetch.
        * \param reduce is an instance of lambda function for reduction.
        * \param keep in an instance of lambda function for storing results.
        * \param zero is zero of given reduction operation also known as idempotent element.
-       * 
+       *
        * \par Example
        * \include Matrices/MultidiagonalMatrix/MultidiagonalMatrixViewExample_rowsReduction.cpp
        * \par Output
@@ -402,7 +410,7 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Method for performing general reduction on all matrix rows for constant instances.
-       * 
+       *
        * \tparam Fetch is a type of lambda function for data fetch declared as
        *          `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`.
        *          The return type of this lambda can be any non void.
@@ -411,12 +419,12 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
        * \tparam Keep is a type of lambda function for storing results of reduction in each row.
        *          It is declared as `keep( const IndexType rowIdx, const double& value )`.
        * \tparam FetchValue is type returned by the Fetch lambda function.
-       * 
+       *
        * \param fetch is an instance of lambda function for data fetch.
        * \param reduce is an instance of lambda function for reduction.
        * \param keep in an instance of lambda function for storing results.
        * \param zero is zero of given reduction operation also known as idempotent element.
-       * 
+       *
        * \par Example
        * \include Matrices/MultidiagonalMatrix/MultidiagonalMatrixViewExample_allRowsReduction.cpp
        * \par Output
@@ -452,7 +460,7 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Method for iteration over all matrix rows for constant instances.
-       * 
+       *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
        *
@@ -486,10 +494,10 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Method for iteration over all matrix rows for non-constant instances.
-       * 
+       *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
-       * 
+       *
        *  `function( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value, bool& compute )`,
        *
        * where
@@ -520,12 +528,12 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief This method calls \e forRows for all matrix rows (for constant instances).
-       * 
+       *
        * See \ref MultidiagonalMatrix::forRows.
-       * 
+       *
        * \tparam Function is a type of lambda function that will operate on matrix elements.
        * \param function  is an instance of the lambda function to be called in each row.
-       * 
+       *
        * \par Example
        * \include Matrices/MultidiagonalMatrix/MultidiagonalMatrixViewExample_forAllRows.cpp
        * \par Output
@@ -536,12 +544,12 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief This method calls \e forRows for all matrix rows.
-       * 
+       *
        * See \ref MultidiagonalMatrix::forRows.
-       * 
+       *
        * \tparam Function is a type of lambda function that will operate on matrix elements.
        * \param function  is an instance of the lambda function to be called in each row.
-       * 
+       *
        * \par Example
        * \include Matrices/MultidiagonalMatrix/MultidiagonalMatrixViewExample_forAllRows.cpp
        * \par Output
@@ -608,16 +616,16 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Computes product of matrix and vector.
-       * 
+       *
        * More precisely, it computes:
-       * 
+       *
        * `outVector = matrixMultiplicator * ( * this ) * inVector + outVectorMultiplicator * outVector`
-       * 
+       *
        * \tparam InVector is type of input vector.  It can be \ref Vector,
        *     \ref VectorView, \ref Array, \ref ArraView or similar container.
        * \tparam OutVector is type of output vector. It can be \ref Vector,
        *     \ref VectorView, \ref Array, \ref ArraView or similar container.
-       * 
+       *
        * \param inVector is input vector.
        * \param outVector is output vector.
        * \param matrixMultiplicator is a factor by which the matrix is multiplied. It is one by default.
@@ -655,7 +663,7 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Assignment of exactly the same matrix type.
-       * 
+       *
        * \param matrix is input matrix for the assignment.
        * \return reference to this matrix.
        */
@@ -663,28 +671,28 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Method for saving the matrix to a file.
-       * 
+       *
        * \param file is the output file.
        */
       void save( File& file ) const;
 
       /**
        * \brief Method for saving the matrix to the file with given filename.
-       * 
+       *
        * \param fileName is name of the file.
        */
       void save( const String& fileName ) const;
 
       /**
        * \brief Method for printing the matrix to output stream.
-       * 
+       *
        * \param str is the output stream.
        */
       void print( std::ostream& str ) const;
 
       /**
        * \brief This method returns matrix elements indexer used by this matrix.
-       * 
+       *
        * \return constant reference to the indexer.
        */
       __cuda_callable__
@@ -692,7 +700,7 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief This method returns matrix elements indexer used by this matrix.
-       * 
+       *
        * \return non-constant reference to the indexer.
        */
       __cuda_callable__
@@ -700,9 +708,9 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Returns padding index denoting padding zero elements.
-       * 
+       *
        * These elements are used for efficient data alignment in memory.
-       * 
+       *
        * \return value of the padding index.
        */
       __cuda_callable__
diff --git a/src/TNL/Matrices/MultidiagonalMatrixView.hpp b/src/TNL/Matrices/MultidiagonalMatrixView.hpp
index 8d1a1d4fe2..1b0687a8cb 100644
--- a/src/TNL/Matrices/MultidiagonalMatrixView.hpp
+++ b/src/TNL/Matrices/MultidiagonalMatrixView.hpp
@@ -101,11 +101,29 @@ template< typename Real,
           typename Index,
           ElementsOrganization Organization >
 __cuda_callable__
-const Index&
+const Index
 MultidiagonalMatrixView< Real, Device, Index, Organization >::
 getDiagonalsCount() const
 {
+#ifdef __CUDA_ARCH__
    return this->diagonalsOffsets.getSize();
+#else
+   return this->hostDiagonalsOffsets.getSize();
+#endif
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          ElementsOrganization Organization >
+   template< typename Vector >
+void
+MultidiagonalMatrixView< Real, Device, Index, Organization >::
+getRowCapacities( Vector& rowCapacities ) const
+{
+   rowCapacities.setSize( this->getRows() );
+   auto aux = this->getDiagonalsCount();
+   rowCapacities = aux;
 }
 
 template< typename Real,
diff --git a/src/TNL/Matrices/SparseMatrix.h b/src/TNL/Matrices/SparseMatrix.h
index 69e02d3c89..247ba8ef42 100644
--- a/src/TNL/Matrices/SparseMatrix.h
+++ b/src/TNL/Matrices/SparseMatrix.h
@@ -398,6 +398,16 @@ class SparseMatrix : public Matrix< Real, Device, Index, RealAllocator >
       template< typename RowsCapacitiesVector >
       void setRowCapacities( const RowsCapacitiesVector& rowCapacities );
 
+      /**
+       * \brief Compute capacities of all rows.
+       *
+       * The row capacities are not stored explicitly and must be computed.
+       *
+       * \param rowCapacities is a vector where the row capacities will be stored.
+       */
+      template< typename Vector >
+      void getRowCapacities( Vector& rowCapacities ) const;
+
       /**
        * \brief This method sets the sparse matrix elements from initializer list.
        *
@@ -452,6 +462,7 @@ class SparseMatrix : public Matrix< Real, Device, Index, RealAllocator >
       template< typename Vector >
       void getCompressedRowLengths( Vector& rowLengths ) const;
 
+
       /**
        * \brief Returns capacity of given matrix row.
        *
diff --git a/src/TNL/Matrices/SparseMatrix.hpp b/src/TNL/Matrices/SparseMatrix.hpp
index 856d52983a..b895fd024a 100644
--- a/src/TNL/Matrices/SparseMatrix.hpp
+++ b/src/TNL/Matrices/SparseMatrix.hpp
@@ -270,6 +270,22 @@ setRowCapacities( const RowsCapacitiesVector& rowsCapacities )
    this->view = this->getView();
 }
 
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
+          typename ComputeReal,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Vector >
+void
+SparseMatrix< Real, Device, Index, MatrixType, Segments, ComputeReal, RealAllocator, IndexAllocator >::
+getRowCapacities( Vector& rowCapacities ) const
+{
+   this->view.getRowCapacities( rowCapacities );
+}
+
 template< typename Real,
           typename Device,
           typename Index,
@@ -905,10 +921,10 @@ operator=( const RHSMatrix& matrix )
    using RHSDeviceType = typename RHSMatrix::DeviceType;
    using RHSRealAllocatorType = typename RHSMatrix::RealAllocatorType;
 
-   Containers::Vector< RHSIndexType, RHSDeviceType, RHSIndexType > rowLengths;
-   matrix.getCompressedRowLengths( rowLengths );
+   Containers::Vector< RHSIndexType, RHSDeviceType, RHSIndexType > rowCapacities;
+   matrix.getRowCapacities( rowCapacities );
    this->setDimensions( matrix.getRows(), matrix.getColumns() );
-   this->setRowCapacities( rowLengths );
+   this->setRowCapacities( rowCapacities );
    Containers::Vector< IndexType, DeviceType, IndexType > rowLocalIndexes( matrix.getRows() );
    rowLocalIndexes = 0;
 
@@ -938,7 +954,7 @@ operator=( const RHSMatrix& matrix )
    }
    else
    {
-      const IndexType maxRowLength = max( rowLengths );
+      const IndexType maxRowLength = max( rowCapacities );
       const IndexType bufferRowsCount( 128 );
       const size_t bufferSize = bufferRowsCount * maxRowLength;
       Containers::Vector< RHSRealType, RHSDeviceType, RHSIndexType, RHSRealAllocatorType > matrixValuesBuffer( bufferSize );
@@ -946,7 +962,9 @@ operator=( const RHSMatrix& matrix )
       Containers::Vector< RealType, DeviceType, IndexType, RealAllocatorType > thisValuesBuffer( bufferSize );
       Containers::Vector< IndexType, DeviceType, IndexType > thisColumnsBuffer( bufferSize );
       Containers::Vector< IndexType, DeviceType, IndexType > thisRowLengths;
-      thisRowLengths = rowLengths;
+      Containers::Vector< RHSIndexType, RHSDeviceType, RHSIndexType > rhsRowLengths;
+      matrix.getCompressedRowLengths( rhsRowLengths );
+      thisRowLengths= rhsRowLengths;
       auto matrixValuesBuffer_view = matrixValuesBuffer.getView();
       auto matrixColumnsBuffer_view = matrixColumnsBuffer.getView();
       auto thisValuesBuffer_view = thisValuesBuffer.getView();
@@ -966,7 +984,11 @@ operator=( const RHSMatrix& matrix )
          auto f1 = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx, RHSIndexType columnIndex, const RHSRealType& value, bool& compute ) mutable {
             if( columnIndex != paddingIndex )
             {
+               //printf("SparseMatrix.hpp: localIdx = %d, maxRowLength = %d \n", localIdx, maxRowLength );
+               TNL_ASSERT_LT( rowIdx - baseRow, bufferRowsCount, "" );
+               TNL_ASSERT_LT( localIdx, maxRowLength, "" );
                const IndexType bufferIdx = ( rowIdx - baseRow ) * maxRowLength + localIdx;
+               TNL_ASSERT_LT( bufferIdx, bufferSize, "" );
                matrixColumnsBuffer_view[ bufferIdx ] = columnIndex;
                matrixValuesBuffer_view[ bufferIdx ] = value;
                //printf( "TO BUFFER: rowIdx = %d localIdx = %d bufferIdx = %d column = %d value = %d \n", rowIdx, localIdx, bufferIdx, columnIndex, value );
diff --git a/src/TNL/Matrices/SparseMatrixView.h b/src/TNL/Matrices/SparseMatrixView.h
index 77dc11f157..24a23f4b6e 100644
--- a/src/TNL/Matrices/SparseMatrixView.h
+++ b/src/TNL/Matrices/SparseMatrixView.h
@@ -249,6 +249,16 @@ class SparseMatrixView : public MatrixView< Real, Device, Index >
       template< typename Vector >
       void getCompressedRowLengths( Vector& rowLengths ) const;
 
+      /**
+       * \brief Compute capacities of all rows.
+       *
+       * The row capacities are not stored explicitly and must be computed.
+       *
+       * \param rowCapacities is a vector where the row capacities will be stored.
+       */
+      template< typename Vector >
+      void getRowCapacities( Vector& rowCapacities ) const;
+
       /**
        * \brief Returns capacity of given matrix row.
        *
@@ -329,7 +339,7 @@ class SparseMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Add element at given \e row and \e column to given \e value.
-       * 
+       *
        * This method can be called from the host system (CPU) no matter
        * where the matrix is allocated. If the matrix is allocated on GPU this method
        * can be called even from device kernels. If the matrix is allocated in GPU device
@@ -337,18 +347,17 @@ class SparseMatrixView : public MatrixView< Real, Device, Index >
        * performance is very low. For higher performance see. \ref SparseMatrix::getRow
        * or \ref SparseMatrix::forRows and \ref SparseMatrix::forAllRows.
        * The call may fail if the matrix row capacity is exhausted.
-       * 
+       *
        * \param row is row index of the element.
        * \param column is columns index of the element.
        * \param value is the value the element will be set to.
        * \param thisElementMultiplicator is multiplicator the original matrix element
        *   value is multiplied by before addition of given \e value.
-       * 
+       *
        * \par Example
        * \include Matrices/SparseMatrix/SparseMatrixViewExample_addElement.cpp
        * \par Output
        * \include SparseMatrixViewExample_addElement.out
-       * 
        */
       __cuda_callable__
       void addElement( IndexType row,
@@ -358,24 +367,24 @@ class SparseMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Returns value of matrix element at position given by its row and column index.
-       * 
+       *
        * This method can be called from the host system (CPU) no matter
        * where the matrix is allocated. If the matrix is allocated on GPU this method
        * can be called even from device kernels. If the matrix is allocated in GPU device
        * this method is called from CPU, it transfers values of each matrix element separately and so the
        * performance is very low. For higher performance see. \ref SparseMatrix::getRow
        * or \ref SparseMatrix::forRows and \ref SparseMatrix::forAllRows.
-       * 
+       *
        * \param row is a row index of the matrix element.
        * \param column i a column index of the matrix element.
-       * 
+       *
        * \return value of given matrix element.
-       * 
+       *
        * \par Example
        * \include Matrices/SparseMatrix/SparseMatrixViewExample_getElement.cpp
        * \par Output
        * \include SparseMatrixViewExample_getElement.out
-       * 
+       *
        */
       __cuda_callable__
       RealType getElement( IndexType row,
@@ -383,7 +392,7 @@ class SparseMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Method for performing general reduction on matrix rows.
-       * 
+       *
        * \tparam Fetch is a type of lambda function for data fetch declared as
        *          `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`.
        *          The return type of this lambda can be any non void.
@@ -392,14 +401,14 @@ class SparseMatrixView : public MatrixView< Real, Device, Index >
        * \tparam Keep is a type of lambda function for storing results of reduction in each row.
        *          It is declared as `keep( const IndexType rowIdx, const double& value )`.
        * \tparam FetchValue is type returned by the Fetch lambda function.
-       * 
+       *
        * \param begin defines beginning of the range [begin,end) of rows to be processed.
        * \param end defines ending of the range [begin,end) of rows to be processed.
        * \param fetch is an instance of lambda function for data fetch.
        * \param reduce is an instance of lambda function for reduction.
        * \param keep in an instance of lambda function for storing results.
        * \param zero is zero of given reduction operation also known as idempotent element.
-       * 
+       *
        * \par Example
        * \include Matrices/SparseMatrix/SparseMatrixViewExample_rowsReduction.cpp
        * \par Output
@@ -410,7 +419,7 @@ class SparseMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Method for performing general reduction on matrix rows for constant instances.
-       * 
+       *
        * \tparam Fetch is a type of lambda function for data fetch declared as
        *          `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`.
        *          The return type of this lambda can be any non void.
@@ -419,14 +428,14 @@ class SparseMatrixView : public MatrixView< Real, Device, Index >
        * \tparam Keep is a type of lambda function for storing results of reduction in each row.
        *          It is declared as `keep( const IndexType rowIdx, const double& value )`.
        * \tparam FetchValue is type returned by the Fetch lambda function.
-       * 
+       *
        * \param begin defines beginning of the range [begin,end) of rows to be processed.
        * \param end defines ending of the range [begin,end) of rows to be processed.
        * \param fetch is an instance of lambda function for data fetch.
        * \param reduce is an instance of lambda function for reduction.
        * \param keep in an instance of lambda function for storing results.
        * \param zero is zero of given reduction operation also known as idempotent element.
-       * 
+       *
        * \par Example
        * \include Matrices/SparseMatrix/SparseMatrixViewExample_rowsReduction.cpp
        * \par Output
@@ -437,7 +446,7 @@ class SparseMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Method for performing general reduction on all matrix rows.
-       * 
+       *
        * \tparam Fetch is a type of lambda function for data fetch declared as
        *          `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`.
        *          The return type of this lambda can be any non void.
@@ -446,12 +455,12 @@ class SparseMatrixView : public MatrixView< Real, Device, Index >
        * \tparam Keep is a type of lambda function for storing results of reduction in each row.
        *          It is declared as `keep( const IndexType rowIdx, const double& value )`.
        * \tparam FetchValue is type returned by the Fetch lambda function.
-       * 
+       *
        * \param fetch is an instance of lambda function for data fetch.
        * \param reduce is an instance of lambda function for reduction.
        * \param keep in an instance of lambda function for storing results.
        * \param zero is zero of given reduction operation also known as idempotent element.
-       * 
+       *
        * \par Example
        * \include Matrices/SparseMatrix/SparseMatrixViewExample_allRowsReduction.cpp
        * \par Output
@@ -462,7 +471,7 @@ class SparseMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Method for performing general reduction on all matrix rows for constant instances.
-       * 
+       *
        * \tparam Fetch is a type of lambda function for data fetch declared as
        *          `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`.
        *          The return type of this lambda can be any non void.
@@ -471,12 +480,12 @@ class SparseMatrixView : public MatrixView< Real, Device, Index >
        * \tparam Keep is a type of lambda function for storing results of reduction in each row.
        *          It is declared as `keep( const IndexType rowIdx, const double& value )`.
        * \tparam FetchValue is type returned by the Fetch lambda function.
-       * 
+       *
        * \param fetch is an instance of lambda function for data fetch.
        * \param reduce is an instance of lambda function for reduction.
        * \param keep in an instance of lambda function for storing results.
        * \param zero is zero of given reduction operation also known as idempotent element.
-       * 
+       *
        * \par Example
        * \include Matrices/SparseMatrix/SparseMatrixViewExample_allRowsReduction.cpp
        * \par Output
@@ -487,18 +496,18 @@ class SparseMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Method for iteration over all matrix rows for constant instances.
-       * 
+       *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
        *  `function( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value, bool& compute )`.
-       *  The \e localIdx parameter is a rank of the non-zero element in given row. 
-       *  If the 'compute' variable is set to false the iteration over the row can 
+       *  The \e localIdx parameter is a rank of the non-zero element in given row.
+       *  If the 'compute' variable is set to false the iteration over the row can
        *  be interrupted.
-       * 
+       *
        * \param begin defines beginning of the range [begin,end) of rows to be processed.
        * \param end defines ending of the range [begin,end) of rows to be processed.
        * \param function is an instance of the lambda function to be called in each row.
-       * 
+       *
        * \par Example
        * \include Matrices/SparseMatrix/SparseMatrixViewExample_forRows.cpp
        * \par Output
@@ -509,18 +518,18 @@ class SparseMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Method for iteration over all matrix rows for non-constant instances.
-       * 
+       *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
        *  `function( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value, bool& compute )`.
-       *  The \e localIdx parameter is a rank of the non-zero element in given row. 
-       *  If the 'compute' variable is set to false the iteration over the row can 
+       *  The \e localIdx parameter is a rank of the non-zero element in given row.
+       *  If the 'compute' variable is set to false the iteration over the row can
        *  be interrupted.
-       * 
+       *
        * \param begin defines beginning of the range [begin,end) of rows to be processed.
        * \param end defines ending of the range [begin,end) of rows to be processed.
        * \param function is an instance of the lambda function to be called in each row.
-       * 
+       *
        * \par Example
        * \include Matrices/SparseMatrix/SparseMatrixViewExample_forRows.cpp
        * \par Output
@@ -531,12 +540,12 @@ class SparseMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief This method calls \e forRows for all matrix rows (for constant instances).
-       * 
+       *
        * See \ref SparseMatrix::forRows.
-       * 
+       *
        * \tparam Function is a type of lambda function that will operate on matrix elements.
        * \param function  is an instance of the lambda function to be called in each row.
-       * 
+       *
        * \par Example
        * \include Matrices/SparseMatrix/SparseMatrixViewExample_forAllRows.cpp
        * \par Output
@@ -547,12 +556,12 @@ class SparseMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief This method calls \e forRows for all matrix rows.
-       * 
+       *
        * See \ref SparseMatrix::forRows.
-       * 
+       *
        * \tparam Function is a type of lambda function that will operate on matrix elements.
        * \param function  is an instance of the lambda function to be called in each row.
-       * 
+       *
        * \par Example
        * \include Matrices/SparseMatrix/SparseMatrixViewExample_forAllRows.cpp
        * \par Output
diff --git a/src/TNL/Matrices/SparseMatrixView.hpp b/src/TNL/Matrices/SparseMatrixView.hpp
index 1a55a22c43..08ac0143a8 100644
--- a/src/TNL/Matrices/SparseMatrixView.hpp
+++ b/src/TNL/Matrices/SparseMatrixView.hpp
@@ -137,6 +137,29 @@ getCompressedRowLengths( Vector& rowLengths ) const
    this->allRowsReduction( fetch, std::plus<>{}, keep, 0 );
 }
 
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          template< typename, typename > class SegmentsView,
+          typename ComputeReal >
+   template< typename Vector >
+void
+SparseMatrixView< Real, Device, Index, MatrixType, SegmentsView, ComputeReal >::
+getRowCapacities( Vector& rowLengths ) const
+{
+   details::set_size_if_resizable( rowLengths, this->getRows() );
+   rowLengths = 0;
+   auto rowLengths_view = rowLengths.getView();
+   auto fetch = [] __cuda_callable__ ( IndexType row, IndexType column, const RealType& value ) -> IndexType {
+      return 1;
+   };
+   auto keep = [=] __cuda_callable__ ( const IndexType rowIdx, const IndexType value ) mutable {
+      rowLengths_view[ rowIdx ] = value;
+   };
+   this->allRowsReduction( fetch, std::plus<>{}, keep, 0 );
+}
+
 template< typename Real,
           typename Device,
           typename Index,
diff --git a/src/TNL/Matrices/TridiagonalMatrix.h b/src/TNL/Matrices/TridiagonalMatrix.h
index 358002ed3b..5a28a34a60 100644
--- a/src/TNL/Matrices/TridiagonalMatrix.h
+++ b/src/TNL/Matrices/TridiagonalMatrix.h
@@ -269,6 +269,16 @@ class TridiagonalMatrix : public Matrix< Real, Device, Index, RealAllocator >
       template< typename ListReal >
       void setElements( const std::initializer_list< std::initializer_list< ListReal > >& data );
 
+      /**
+       * \brief Compute capacities of all rows.
+       *
+       * The row capacities are not stored explicitly and must be computed.
+       *
+       * \param rowCapacities is a vector where the row capacities will be stored.
+       */
+      template< typename Vector >
+      void getRowCapacities( Vector& rowCapacities ) const;
+
       /**
        * \brief Computes number of non-zeros in each row.
        *
diff --git a/src/TNL/Matrices/TridiagonalMatrix.hpp b/src/TNL/Matrices/TridiagonalMatrix.hpp
index 37d1f14507..a6f511470e 100644
--- a/src/TNL/Matrices/TridiagonalMatrix.hpp
+++ b/src/TNL/Matrices/TridiagonalMatrix.hpp
@@ -135,6 +135,19 @@ setRowCapacities( const RowCapacitiesVector& rowCapacities )
          throw std::logic_error( "Too many non-zero elements per row in a tri-diagonal matrix." );
 }
 
+template< typename Real,
+          typename Device,
+          typename Index,
+          ElementsOrganization Organization,
+          typename RealAllocator >
+   template< typename Vector >
+void
+TridiagonalMatrix< Real, Device, Index, Organization, RealAllocator >::
+getRowCapacities( Vector& rowCapacities ) const
+{
+   return this->view.getRowCapacities( rowCapacities );
+}
+
 template< typename Real,
           typename Device,
           typename Index,
diff --git a/src/TNL/Matrices/TridiagonalMatrixView.h b/src/TNL/Matrices/TridiagonalMatrixView.h
index 4fc6c86cd2..49eeec3b6d 100644
--- a/src/TNL/Matrices/TridiagonalMatrixView.h
+++ b/src/TNL/Matrices/TridiagonalMatrixView.h
@@ -63,7 +63,6 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
       using IndexType = Index;
 
       /**
-       * \brief Type of related matrix view. 
        */
       using ViewType = TridiagonalMatrixView< Real, Device, Index, Organization >;
 
@@ -94,7 +93,7 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Constructor with all necessary data and views.
-       * 
+       *
        * \param values is a vector view with matrix elements values
        * \param indexer is an indexer of matrix elements
        */
@@ -103,7 +102,7 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Copy constructor.
-       * 
+       *
        * \param matrix is an input tridiagonal matrix view.
        */
       __cuda_callable__
@@ -111,7 +110,7 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Move constructor.
-       * 
+       *
        * \param matrix is an input tridiagonal matrix view.
        */
       __cuda_callable__
@@ -119,44 +118,54 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Returns a modifiable view of the tridiagonal matrix.
-       * 
+       *
        * \return tridiagonal matrix view.
        */
       ViewType getView();
 
       /**
        * \brief Returns a non-modifiable view of the tridiagonal matrix.
-       * 
+       *
        * \return tridiagonal matrix view.
        */
       ConstViewType getConstView() const;
 
       /**
        * \brief Returns string with serialization type.
-       * 
+       *
        * The string has a form `Matrices::TridiagonalMatrix< RealType,  [any_device], IndexType, Organization, [any_allocator] >`.
-       * 
+       *
        * See \ref TridiagonalMatrix::getSerializationType.
-       * 
+       *
        * \return \ref String with the serialization type.
        */
       static String getSerializationType();
 
       /**
        * \brief Returns string with serialization type.
-       * 
+       *
        * See \ref TridiagonalMatrix::getSerializationType.
-       * 
+       *
        * \return \ref String with the serialization type.
        */
       virtual String getSerializationTypeVirtual() const;
 
+      /**
+       * \brief Compute capacities of all rows.
+       *
+       * The row capacities are not stored explicitly and must be computed.
+       *
+       * \param rowCapacities is a vector where the row capacities will be stored.
+       */
+      template< typename Vector >
+      void getRowCapacities( Vector& rowCapacities ) const;
+
       /**
        * \brief Computes number of non-zeros in each row.
-       * 
+       *
        * \param rowLengths is a vector into which the number of non-zeros in each row
        * will be stored.
-       * 
+       *
        * \par Example
        * \include Matrices/TridiagonalMatrix/TridiagonalMatrixViewExample_getCompressedRowLengths.cpp
        * \par Output
@@ -182,12 +191,12 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Comparison operator with another tridiagonal matrix.
-       * 
+       *
        * \tparam Real_ is \e Real type of the source matrix.
        * \tparam Device_ is \e Device type of the source matrix.
        * \tparam Index_ is \e Index type of the source matrix.
        * \tparam Organization_ is \e Organization of the source matrix.
-       * 
+       *
        * \return \e true if both matrices are identical and \e false otherwise.
        */
       template< typename Real_,
@@ -198,14 +207,14 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Comparison operator with another multidiagonal matrix.
-       * 
+       *
        * \tparam Real_ is \e Real type of the source matrix.
        * \tparam Device_ is \e Device type of the source matrix.
        * \tparam Index_ is \e Index type of the source matrix.
        * \tparam Organization_ is \e Organization of the source matrix.
-       * 
+       *
        * \param matrix is the source matrix.
-       * 
+       *
        * \return \e true if both matrices are NOT identical and \e false otherwise.
        */
       template< typename Real_,
@@ -216,16 +225,16 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Non-constant getter of simple structure for accessing given matrix row.
-       * 
+       *
        * \param rowIdx is matrix row index.
-       * 
+       *
        * \return RowView for accessing given matrix row.
-       * 
+       *
        * \par Example
        * \include Matrices/TridiagonalMatrix/TridiagonalMatrixViewExample_getRow.cpp
        * \par Output
        * \include TridiagonalMatrixViewExample_getRow.out
-       * 
+       *
        * See \ref TridiagonalMatrixRowView.
        */
       __cuda_callable__
@@ -233,16 +242,16 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Constant getter of simple structure for accessing given matrix row.
-       * 
+       *
        * \param rowIdx is matrix row index.
-       * 
+       *
        * \return RowView for accessing given matrix row.
        *
        * \par Example
        * \include Matrices/TridiagonalMatrix/TridiagonalMatrixViewExample_getConstRow.cpp
        * \par Output
        * \include TridiagonalMatrixViewExample_getConstRow.out
-       * 
+       *
        * See \ref TridiagonalMatrixRowView.
        */
       __cuda_callable__
@@ -250,14 +259,14 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Set all matrix elements to given value.
-       * 
+       *
        * \param value is the new value of all matrix elements.
        */
       void setValue( const RealType& v );
 
       /**
        * \brief Sets element at given \e row and \e column to given \e value.
-       * 
+       *
        * This method can be called from the host system (CPU) no matter
        * where the matrix is allocated. If the matrix is allocated on GPU this method
        * can be called even from device kernels. If the matrix is allocated in GPU device
@@ -265,11 +274,11 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
        * performance is very low. For higher performance see. \ref TridiagonalMatrix::getRow
        * or \ref TridiagonalMatrix::forRows and \ref TridiagonalMatrix::forAllRows.
        * The call may fail if the matrix row capacity is exhausted.
-       * 
+       *
        * \param row is row index of the element.
        * \param column is columns index of the element.
        * \param value is the value the element will be set to.
-       * 
+       *
        * \par Example
        * \include Matrices/TridiagonalMatrix/TridiagonalMatrixViewExample_setElement.cpp
        * \par Output
@@ -282,7 +291,7 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Add element at given \e row and \e column to given \e value.
-       * 
+       *
        * This method can be called from the host system (CPU) no matter
        * where the matrix is allocated. If the matrix is allocated on GPU this method
        * can be called even from device kernels. If the matrix is allocated in GPU device
@@ -290,18 +299,17 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
        * performance is very low. For higher performance see. \ref TridiagonalMatrix::getRow
        * or \ref TridiagonalMatrix::forRows and \ref TridiagonalMatrix::forAllRows.
        * The call may fail if the matrix row capacity is exhausted.
-       * 
+       *
        * \param row is row index of the element.
        * \param column is columns index of the element.
        * \param value is the value the element will be set to.
        * \param thisElementMultiplicator is multiplicator the original matrix element
        *   value is multiplied by before addition of given \e value.
-       * 
+       *
        * \par Example
        * \include Matrices/TridiagonalMatrix/TridiagonalMatrixViewExample_addElement.cpp
        * \par Output
        * \include TridiagonalMatrixViewExample_addElement.out
-       * 
        */
       __cuda_callable__
       void addElement( const IndexType row,
@@ -311,24 +319,23 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Returns value of matrix element at position given by its row and column index.
-       * 
+       *
        * This method can be called from the host system (CPU) no matter
        * where the matrix is allocated. If the matrix is allocated on GPU this method
        * can be called even from device kernels. If the matrix is allocated in GPU device
        * this method is called from CPU, it transfers values of each matrix element separately and so the
        * performance is very low. For higher performance see. \ref TridiagonalMatrix::getRow
        * or \ref TridiagonalMatrix::forRows and \ref TridiagonalMatrix::forAllRows.
-       * 
+       *
        * \param row is a row index of the matrix element.
        * \param column i a column index of the matrix element.
-       * 
+       *
        * \return value of given matrix element.
-       * 
+       *
        * \par Example
        * \include Matrices/TridiagonalMatrix/TridiagonalMatrixViewExample_getElement.cpp
        * \par Output
        * \include TridiagonalMatrixViewExample_getElement.out
-       * 
        */
       __cuda_callable__
       RealType getElement( const IndexType row,
@@ -336,7 +343,7 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Method for performing general reduction on matrix rows for constant instances.
-       * 
+       *
        * \tparam Fetch is a type of lambda function for data fetch declared as
        *          `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`.
        *          The return type of this lambda can be any non void.
@@ -345,14 +352,14 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
        * \tparam Keep is a type of lambda function for storing results of reduction in each row.
        *          It is declared as `keep( const IndexType rowIdx, const double& value )`.
        * \tparam FetchValue is type returned by the Fetch lambda function.
-       * 
+       *
        * \param begin defines beginning of the range [begin,end) of rows to be processed.
        * \param end defines ending of the range [begin,end) of rows to be processed.
        * \param fetch is an instance of lambda function for data fetch.
        * \param reduce is an instance of lambda function for reduction.
        * \param keep in an instance of lambda function for storing results.
        * \param zero is zero of given reduction operation also known as idempotent element.
-       * 
+       *
        * \par Example
        * \include Matrices/TridiagonalMatrix/TridiagonalMatrixViewExample_rowsReduction.cpp
        * \par Output
@@ -363,7 +370,7 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Method for performing general reduction on matrix rows.
-       * 
+       *
        * \tparam Fetch is a type of lambda function for data fetch declared as
        *          `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`.
        *          The return type of this lambda can be any non void.
@@ -372,14 +379,14 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
        * \tparam Keep is a type of lambda function for storing results of reduction in each row.
        *          It is declared as `keep( const IndexType rowIdx, const double& value )`.
        * \tparam FetchValue is type returned by the Fetch lambda function.
-       * 
+       *
        * \param begin defines beginning of the range [begin,end) of rows to be processed.
        * \param end defines ending of the range [begin,end) of rows to be processed.
        * \param fetch is an instance of lambda function for data fetch.
        * \param reduce is an instance of lambda function for reduction.
        * \param keep in an instance of lambda function for storing results.
        * \param zero is zero of given reduction operation also known as idempotent element.
-       * 
+       *
        * \par Example
        * \include Matrices/TridiagonalMatrix/TridiagonalMatrixViewExample_rowsReduction.cpp
        * \par Output
@@ -390,7 +397,7 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Method for performing general reduction on all matrix rows for constant instances.
-       * 
+       *
        * \tparam Fetch is a type of lambda function for data fetch declared as
        *          `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`.
        *          The return type of this lambda can be any non void.
@@ -399,12 +406,12 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
        * \tparam Keep is a type of lambda function for storing results of reduction in each row.
        *          It is declared as `keep( const IndexType rowIdx, const double& value )`.
        * \tparam FetchValue is type returned by the Fetch lambda function.
-       * 
+       *
        * \param fetch is an instance of lambda function for data fetch.
        * \param reduce is an instance of lambda function for reduction.
        * \param keep in an instance of lambda function for storing results.
        * \param zero is zero of given reduction operation also known as idempotent element.
-       * 
+       *
        * \par Example
        * \include Matrices/TridiagonalMatrix/TridiagonalMatrixViewExample_allRowsReduction.cpp
        * \par Output
@@ -415,7 +422,7 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Method for performing general reduction on all matrix rows.
-       * 
+       *
        * \tparam Fetch is a type of lambda function for data fetch declared as
        *          `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`.
        *          The return type of this lambda can be any non void.
@@ -424,12 +431,12 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
        * \tparam Keep is a type of lambda function for storing results of reduction in each row.
        *          It is declared as `keep( const IndexType rowIdx, const double& value )`.
        * \tparam FetchValue is type returned by the Fetch lambda function.
-       * 
+       *
        * \param fetch is an instance of lambda function for data fetch.
        * \param reduce is an instance of lambda function for reduction.
        * \param keep in an instance of lambda function for storing results.
        * \param zero is zero of given reduction operation also known as idempotent element.
-       * 
+       *
        * \par Example
        * \include Matrices/TridiagonalMatrix/TridiagonalMatrixViewExample_allRowsReduction.cpp
        * \par Output
@@ -440,18 +447,18 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Method for iteration over all matrix rows for constant instances.
-       * 
+       *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
        *  `function( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value, bool& compute )`.
-       *  The \e localIdx parameter is a rank of the non-zero element in given row. 
-       *  If the 'compute' variable is set to false the iteration over the row can 
+       *  The \e localIdx parameter is a rank of the non-zero element in given row.
+       *  If the 'compute' variable is set to false the iteration over the row can
        *  be interrupted.
-       * 
+       *
        * \param begin defines beginning of the range [begin,end) of rows to be processed.
        * \param end defines ending of the range [begin,end) of rows to be processed.
        * \param function is an instance of the lambda function to be called in each row.
-       * 
+       *
        * \par Example
        * \include Matrices/TridiagonalMatrix/TridiagonalMatrixViewExample_forRows.cpp
        * \par Output
@@ -462,18 +469,18 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Method for iteration over all matrix rows for non-constant instances.
-       * 
+       *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
        *  `function( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value, bool& compute )`.
-       *  The \e localIdx parameter is a rank of the non-zero element in given row. 
-       *  If the 'compute' variable is set to false the iteration over the row can 
+       *  The \e localIdx parameter is a rank of the non-zero element in given row.
+       *  If the 'compute' variable is set to false the iteration over the row can
        *  be interrupted.
-       * 
+       *
        * \param begin defines beginning of the range [begin,end) of rows to be processed.
        * \param end defines ending of the range [begin,end) of rows to be processed.
        * \param function is an instance of the lambda function to be called in each row.
-       * 
+       *
        * \par Example
        * \include Matrices/TridiagonalMatrix/TridiagonalMatrixViewExample_forRows.cpp
        * \par Output
@@ -484,12 +491,12 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief This method calls \e forRows for all matrix rows (for constant instances).
-       * 
+       *
        * See \ref TridiagonalMatrix::forRows.
-       * 
+       *
        * \tparam Function is a type of lambda function that will operate on matrix elements.
        * \param function  is an instance of the lambda function to be called in each row.
-       * 
+       *
        * \par Example
        * \include Matrices/TridiagonalMatrix/TridiagonalMatrixViewExample_forAllRows.cpp
        * \par Output
@@ -500,12 +507,12 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief This method calls \e forRows for all matrix rows.
-       * 
+       *
        * See \ref TridiagonalMatrix::forRows.
-       * 
+       *
        * \tparam Function is a type of lambda function that will operate on matrix elements.
        * \param function  is an instance of the lambda function to be called in each row.
-       * 
+       *
        * \par Example
        * \include Matrices/TridiagonalMatrix/TridiagonalMatrixViewExample_forAllRows.cpp
        * \par Output
@@ -572,16 +579,16 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Computes product of matrix and vector.
-       * 
+       *
        * More precisely, it computes:
-       * 
+       *
        * `outVector = matrixMultiplicator * ( * this ) * inVector + outVectorMultiplicator * outVector`
-       * 
+       *
        * \tparam InVector is type of input vector.  It can be \ref Vector,
        *     \ref VectorView, \ref Array, \ref ArraView or similar container.
        * \tparam OutVector is type of output vector. It can be \ref Vector,
        *     \ref VectorView, \ref Array, \ref ArraView or similar container.
-       * 
+       *
        * \param inVector is input vector.
        * \param outVector is output vector.
        * \param matrixMultiplicator is a factor by which the matrix is multiplied. It is one by default.
@@ -619,7 +626,7 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Assignment of exactly the same matrix type.
-       * 
+       *
        * \param matrix is input matrix for the assignment.
        * \return reference to this matrix.
        */
@@ -627,28 +634,28 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Method for saving the matrix to a file.
-       * 
+       *
        * \param file is the output file.
        */
       void save( File& file ) const;
 
       /**
        * \brief Method for saving the matrix to the file with given filename.
-       * 
+       *
        * \param fileName is name of the file.
        */
       void save( const String& fileName ) const;
 
       /**
        * \brief Method for printing the matrix to output stream.
-       * 
+       *
        * \param str is the output stream.
        */
       void print( std::ostream& str ) const;
 
       /**
        * \brief This method returns matrix elements indexer used by this matrix.
-       * 
+       *
        * \return constant reference to the indexer.
        */
       __cuda_callable__
@@ -656,7 +663,7 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief This method returns matrix elements indexer used by this matrix.
-       * 
+       *
        * \return non-constant reference to the indexer.
        */
       __cuda_callable__
@@ -664,9 +671,9 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       /**
        * \brief Returns padding index denoting padding zero elements.
-       * 
+       *
        * These elements are used for efficient data alignment in memory.
-       * 
+       *
        * \return value of the padding index.
        */
       __cuda_callable__
diff --git a/src/TNL/Matrices/TridiagonalMatrixView.hpp b/src/TNL/Matrices/TridiagonalMatrixView.hpp
index 0d6bfe064c..595a058cc7 100644
--- a/src/TNL/Matrices/TridiagonalMatrixView.hpp
+++ b/src/TNL/Matrices/TridiagonalMatrixView.hpp
@@ -84,6 +84,19 @@ getSerializationTypeVirtual() const
    return this->getSerializationType();
 }
 
+template< typename Real,
+          typename Device,
+          typename Index,
+          ElementsOrganization Organization >
+   template< typename Vector >
+void
+TridiagonalMatrixView< Real, Device, Index, Organization >::
+getRowCapacities( Vector& rowCapacities ) const
+{
+   rowCapacities.setSize( this->getRows() );
+   rowCapacities = 3;
+}
+
 template< typename Real,
           typename Device,
           typename Index,
-- 
GitLab