From 219abe6c3cc9bce48bdb2391f062f38461bc60cc Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= <oberhuber.tomas@gmail.com>
Date: Tue, 6 Apr 2021 21:44:29 +0200
Subject: [PATCH] Removing boolean variable compute from forElements.

---
 .../SegmentsExample_CSR_constructor_1.cpp     |  2 +-
 .../SegmentsExample_CSR_constructor_2.cpp     |  2 +-
 .../Segments/SegmentsExample_General.cpp      |  2 +-
 .../DenseMatrixExample_forAllElements.cpp     |  6 +--
 .../DenseMatrixExample_forElements.cpp        |  6 +--
 .../DenseMatrixViewExample_forAllElements.cpp |  6 +--
 .../DenseMatrixViewExample_forElements.cpp    |  6 +--
 .../LambdaMatrixExample_forAllElements.cpp    |  2 +-
 .../LambdaMatrixExample_forElements.cpp       |  2 +-
 ...tidiagonalMatrixExample_forAllElements.cpp |  2 +-
 ...MultidiagonalMatrixExample_forElements.cpp |  2 +-
 ...agonalMatrixViewExample_forAllElements.cpp |  2 +-
 ...idiagonalMatrixViewExample_forElements.cpp |  2 +-
 .../SparseMatrixExample_forAllElements.cpp    |  6 +--
 .../SparseMatrixExample_forElements.cpp       |  6 +--
 ...SparseMatrixViewExample_forAllElements.cpp |  6 +--
 .../SparseMatrixViewExample_forElements.cpp   |  6 +--
 ...ridiagonalMatrixExample_forAllElements.cpp |  2 +-
 .../TridiagonalMatrixExample_forElements.cpp  |  2 +-
 ...agonalMatrixViewExample_forAllElements.cpp |  2 +-
 ...idiagonalMatrixViewExample_forElements.cpp |  2 +-
 .../Matrices/DenseMatrixSetup_Benchmark.cpp   |  2 +-
 .../MultidiagonalMatrixSetup_Benchmark.cpp    |  2 +-
 .../Matrices/SparseMatrixSetup_Benchmark.cpp  |  2 +-
 src/TNL/Algorithms/Segments/BiEllpackView.h   | 18 +++------
 src/TNL/Algorithms/Segments/BiEllpackView.hpp |  7 ++--
 src/TNL/Algorithms/Segments/CSR.h             |  2 +-
 src/TNL/Algorithms/Segments/CSRView.hpp       |  5 +--
 .../Algorithms/Segments/ChunkedEllpackView.h  | 18 +++------
 .../Segments/ChunkedEllpackView.hpp           |  9 ++---
 src/TNL/Algorithms/Segments/EllpackView.hpp   | 10 ++---
 .../Algorithms/Segments/SlicedEllpackView.hpp | 14 +++----
 src/TNL/Matrices/DenseMatrix.h                | 16 ++------
 src/TNL/Matrices/DenseMatrix.hpp              |  8 ++--
 src/TNL/Matrices/DenseMatrixView.h            |  8 ++--
 src/TNL/Matrices/DenseMatrixView.hpp          |  8 ++--
 src/TNL/Matrices/LambdaMatrix.h               |  6 +--
 src/TNL/Matrices/LambdaMatrix.hpp             |  5 +--
 src/TNL/Matrices/MultidiagonalMatrix.h        | 18 ++-------
 src/TNL/Matrices/MultidiagonalMatrix.hpp      |  6 +--
 src/TNL/Matrices/MultidiagonalMatrixView.h    | 18 ++-------
 src/TNL/Matrices/MultidiagonalMatrixView.hpp  | 10 ++---
 src/TNL/Matrices/SparseMatrix.h               | 16 ++------
 src/TNL/Matrices/SparseMatrix.hpp             | 12 +++---
 src/TNL/Matrices/SparseMatrixView.h           | 16 ++------
 src/TNL/Matrices/SparseMatrixView.hpp         | 12 +++---
 src/TNL/Matrices/TridiagonalMatrix.h          | 24 +++--------
 src/TNL/Matrices/TridiagonalMatrix.hpp        |  4 +-
 src/TNL/Matrices/TridiagonalMatrixView.h      | 12 ++----
 src/TNL/Matrices/TridiagonalMatrixView.hpp    | 40 +++++++++----------
 src/UnitTests/Matrices/DenseMatrixTest.h      |  2 +-
 src/UnitTests/Matrices/SparseMatrixTest.hpp   |  2 +-
 .../SparseMatrixVectorProductTest.hpp         |  6 +--
 53 files changed, 155 insertions(+), 257 deletions(-)

diff --git a/Documentation/Examples/Algorithms/Segments/SegmentsExample_CSR_constructor_1.cpp b/Documentation/Examples/Algorithms/Segments/SegmentsExample_CSR_constructor_1.cpp
index 4c6e285752..0ceb7a6bd4 100644
--- a/Documentation/Examples/Algorithms/Segments/SegmentsExample_CSR_constructor_1.cpp
+++ b/Documentation/Examples/Algorithms/Segments/SegmentsExample_CSR_constructor_1.cpp
@@ -26,7 +26,7 @@ void SegmentsExample()
     * Insert data into particular segments.
     */
    auto data_view = data.getView();
-   segments.forAllElements( [=] __cuda_callable__ ( int segmentIdx, int localIdx, int globalIdx, bool& compute ) mutable {
+   segments.forAllElements( [=] __cuda_callable__ ( int segmentIdx, int localIdx, int globalIdx ) mutable {
       if( localIdx <= segmentIdx )
          data_view[ globalIdx ] = segmentIdx;
    } );
diff --git a/Documentation/Examples/Algorithms/Segments/SegmentsExample_CSR_constructor_2.cpp b/Documentation/Examples/Algorithms/Segments/SegmentsExample_CSR_constructor_2.cpp
index c15f5791ee..9493758b49 100644
--- a/Documentation/Examples/Algorithms/Segments/SegmentsExample_CSR_constructor_2.cpp
+++ b/Documentation/Examples/Algorithms/Segments/SegmentsExample_CSR_constructor_2.cpp
@@ -25,7 +25,7 @@ void SegmentsExample()
     * Insert data into particular segments.
     */
    auto data_view = data.getView();
-   segments.forAllElements( [=] __cuda_callable__ ( int segmentIdx, int localIdx, int globalIdx, bool& compute ) mutable {
+   segments.forAllElements( [=] __cuda_callable__ ( int segmentIdx, int localIdx, int globalIdx ) mutable {
       if( localIdx <= segmentIdx )
          data_view[ globalIdx ] = segmentIdx;
    } );
diff --git a/Documentation/Examples/Algorithms/Segments/SegmentsExample_General.cpp b/Documentation/Examples/Algorithms/Segments/SegmentsExample_General.cpp
index e50c6d1edd..ade0263fbc 100644
--- a/Documentation/Examples/Algorithms/Segments/SegmentsExample_General.cpp
+++ b/Documentation/Examples/Algorithms/Segments/SegmentsExample_General.cpp
@@ -27,7 +27,7 @@ void SegmentsExample()
     * Insert data into particular segments.
     */
    auto data_view = data.getView();
-   segments.forAllElements( [=] __cuda_callable__ ( IndexType segmentIdx, IndexType localIdx, IndexType globalIdx, bool& compute ) mutable {
+   segments.forAllElements( [=] __cuda_callable__ ( IndexType segmentIdx, IndexType localIdx, IndexType globalIdx ) mutable {
       if( localIdx <= segmentIdx )
          data_view[ globalIdx ] = segmentIdx;
    } );
diff --git a/Documentation/Examples/Matrices/DenseMatrix/DenseMatrixExample_forAllElements.cpp b/Documentation/Examples/Matrices/DenseMatrix/DenseMatrixExample_forAllElements.cpp
index 4fd7d3b472..f143164a25 100644
--- a/Documentation/Examples/Matrices/DenseMatrix/DenseMatrixExample_forAllElements.cpp
+++ b/Documentation/Examples/Matrices/DenseMatrix/DenseMatrixExample_forAllElements.cpp
@@ -8,10 +8,8 @@ void forAllElementsExample()
 {
    TNL::Matrices::DenseMatrix< double, Device > matrix( 5, 5 );
 
-   auto f = [=] __cuda_callable__ ( int rowIdx, int columnIdx, int columnIdx_, double& value, bool& compute ) {
-      if( rowIdx < columnIdx )
-         compute = false;
-      else
+   auto f = [=] __cuda_callable__ ( int rowIdx, int columnIdx, int columnIdx_, double& value ) {
+      if( rowIdx >= columnIdx )
          value = rowIdx + columnIdx;
    };
 
diff --git a/Documentation/Examples/Matrices/DenseMatrix/DenseMatrixExample_forElements.cpp b/Documentation/Examples/Matrices/DenseMatrix/DenseMatrixExample_forElements.cpp
index 0764eecdfc..b37470c43c 100644
--- a/Documentation/Examples/Matrices/DenseMatrix/DenseMatrixExample_forElements.cpp
+++ b/Documentation/Examples/Matrices/DenseMatrix/DenseMatrixExample_forElements.cpp
@@ -8,10 +8,8 @@ void forElementsExample()
 {
    TNL::Matrices::DenseMatrix< double, Device > matrix( 5, 5 );
 
-   auto f = [=] __cuda_callable__ ( int rowIdx, int columnIdx, int columnIdx_, double& value, bool& compute ) {
-      if( rowIdx < columnIdx )
-         compute = false;
-      else
+   auto f = [=] __cuda_callable__ ( int rowIdx, int columnIdx, int columnIdx_, double& value ) {
+      if( rowIdx >= columnIdx )
          value = rowIdx + columnIdx;
    };
 
diff --git a/Documentation/Examples/Matrices/DenseMatrix/DenseMatrixViewExample_forAllElements.cpp b/Documentation/Examples/Matrices/DenseMatrix/DenseMatrixViewExample_forAllElements.cpp
index 66b3941301..c5802a0e1f 100644
--- a/Documentation/Examples/Matrices/DenseMatrix/DenseMatrixViewExample_forAllElements.cpp
+++ b/Documentation/Examples/Matrices/DenseMatrix/DenseMatrixViewExample_forAllElements.cpp
@@ -9,10 +9,8 @@ void forAllElementsExample()
    TNL::Matrices::DenseMatrix< double, Device > matrix( 5, 5 );
    auto matrixView = matrix.getView();
 
-   auto f = [=] __cuda_callable__ ( int rowIdx, int columnIdx, int globalIdx, double& value, bool& compute ) {
-      if( rowIdx < columnIdx )
-         compute = false;
-      else
+   auto f = [=] __cuda_callable__ ( int rowIdx, int columnIdx, int globalIdx, double& value ) {
+      if( rowIdx >= columnIdx )
          value = rowIdx + columnIdx;
    };
 
diff --git a/Documentation/Examples/Matrices/DenseMatrix/DenseMatrixViewExample_forElements.cpp b/Documentation/Examples/Matrices/DenseMatrix/DenseMatrixViewExample_forElements.cpp
index 6a980d23c1..572c526f02 100644
--- a/Documentation/Examples/Matrices/DenseMatrix/DenseMatrixViewExample_forElements.cpp
+++ b/Documentation/Examples/Matrices/DenseMatrix/DenseMatrixViewExample_forElements.cpp
@@ -9,10 +9,8 @@ void forElementsExample()
    TNL::Matrices::DenseMatrix< double, Device > matrix( 5, 5 );
    auto matrixView = matrix.getView();
 
-   auto f = [=] __cuda_callable__ ( int rowIdx, int columnIdx, int globalIdx, double& value, bool& compute ) {
-      if( columnIdx > rowIdx )
-         compute = false;
-      else
+   auto f = [=] __cuda_callable__ ( int rowIdx, int columnIdx, int globalIdx, double& value ) {
+      if( columnIdx <= rowIdx )
          value = rowIdx + columnIdx;
    };
 
diff --git a/Documentation/Examples/Matrices/LambdaMatrix/LambdaMatrixExample_forAllElements.cpp b/Documentation/Examples/Matrices/LambdaMatrix/LambdaMatrixExample_forAllElements.cpp
index 293f173d2a..6b335f5f29 100644
--- a/Documentation/Examples/Matrices/LambdaMatrix/LambdaMatrixExample_forAllElements.cpp
+++ b/Documentation/Examples/Matrices/LambdaMatrix/LambdaMatrixExample_forAllElements.cpp
@@ -22,7 +22,7 @@ void forAllElementsExample()
    TNL::Matrices::DenseMatrix< double, Device > denseMatrix( 5, 5 );
    auto denseView = denseMatrix.getView();
 
-   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int columnIdx, double value, bool& compute ) mutable {
+   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int columnIdx, double value ) mutable {
       denseView.setElement( rowIdx, columnIdx, value );
    };
 
diff --git a/Documentation/Examples/Matrices/LambdaMatrix/LambdaMatrixExample_forElements.cpp b/Documentation/Examples/Matrices/LambdaMatrix/LambdaMatrixExample_forElements.cpp
index f23f031b1b..8472ef28d2 100644
--- a/Documentation/Examples/Matrices/LambdaMatrix/LambdaMatrixExample_forElements.cpp
+++ b/Documentation/Examples/Matrices/LambdaMatrix/LambdaMatrixExample_forElements.cpp
@@ -22,7 +22,7 @@ void forElementsExample()
    TNL::Matrices::DenseMatrix< double, Device > denseMatrix( 5, 5 );
    auto denseView = denseMatrix.getView();
 
-   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int columnIdx, double value, bool& compute ) mutable {
+   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int columnIdx, double value ) mutable {
       denseView.setElement( rowIdx, columnIdx, value );
    };
 
diff --git a/Documentation/Examples/Matrices/MultidiagonalMatrix/MultidiagonalMatrixExample_forAllElements.cpp b/Documentation/Examples/Matrices/MultidiagonalMatrix/MultidiagonalMatrixExample_forAllElements.cpp
index b29543d9e0..96a4668b96 100644
--- a/Documentation/Examples/Matrices/MultidiagonalMatrix/MultidiagonalMatrixExample_forAllElements.cpp
+++ b/Documentation/Examples/Matrices/MultidiagonalMatrix/MultidiagonalMatrixExample_forAllElements.cpp
@@ -23,7 +23,7 @@ void forAllElementsExample()
       5,               // number of matrix columns
       { -2, -1, 0 } ); // matrix diagonals offsets
 
-   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int columnIdx, double& value, bool& compute ) {
+   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int columnIdx, double& value ) {
       /***
        * 'forElements' method iterates only over matrix elements lying on given subdiagonals
        * and so we do not need to check anything. The element value can be expressed
diff --git a/Documentation/Examples/Matrices/MultidiagonalMatrix/MultidiagonalMatrixExample_forElements.cpp b/Documentation/Examples/Matrices/MultidiagonalMatrix/MultidiagonalMatrixExample_forElements.cpp
index dd30694e6f..1dc957af24 100644
--- a/Documentation/Examples/Matrices/MultidiagonalMatrix/MultidiagonalMatrixExample_forElements.cpp
+++ b/Documentation/Examples/Matrices/MultidiagonalMatrix/MultidiagonalMatrixExample_forElements.cpp
@@ -23,7 +23,7 @@ void forElementsExample()
       5,               // number of matrix columns
       { -2, -1, 0 } ); // matrix diagonals offsets
 
-   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int columnIdx, double& value, bool& compute ) {
+   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int columnIdx, double& value ) {
       /***
        * 'forElements' method iterates only over matrix elements lying on given subdiagonals
        * and so we do not need to check anything. The element value can be expressed
diff --git a/Documentation/Examples/Matrices/MultidiagonalMatrix/MultidiagonalMatrixViewExample_forAllElements.cpp b/Documentation/Examples/Matrices/MultidiagonalMatrix/MultidiagonalMatrixViewExample_forAllElements.cpp
index b05da1d821..4ca0940cbe 100644
--- a/Documentation/Examples/Matrices/MultidiagonalMatrix/MultidiagonalMatrixViewExample_forAllElements.cpp
+++ b/Documentation/Examples/Matrices/MultidiagonalMatrix/MultidiagonalMatrixViewExample_forAllElements.cpp
@@ -24,7 +24,7 @@ void forAllElementsExample()
       { -2, -1, 0 } ); // matrix diagonals offsets
    auto view = matrix.getView();
 
-   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int columnIdx, double& value, bool& compute ) {
+   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int columnIdx, double& value ) {
       /***
        * 'forElements' method iterates only over matrix elements lying on given subdiagonals
        * and so we do not need to check anything. The element value can be expressed
diff --git a/Documentation/Examples/Matrices/MultidiagonalMatrix/MultidiagonalMatrixViewExample_forElements.cpp b/Documentation/Examples/Matrices/MultidiagonalMatrix/MultidiagonalMatrixViewExample_forElements.cpp
index 9663a2c0d4..d941fc4a20 100644
--- a/Documentation/Examples/Matrices/MultidiagonalMatrix/MultidiagonalMatrixViewExample_forElements.cpp
+++ b/Documentation/Examples/Matrices/MultidiagonalMatrix/MultidiagonalMatrixViewExample_forElements.cpp
@@ -24,7 +24,7 @@ void forElementsExample()
       { -2, -1, 0 } ); // matrix diagonals offsets
    auto view = matrix.getView();
 
-   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int columnIdx, double& value, bool& compute ) {
+   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int columnIdx, double& value ) {
       /***
        * 'forElements' method iterates only over matrix elements lying on given subdiagonals
        * and so we do not need to check anything. The element value can be expressed
diff --git a/Documentation/Examples/Matrices/SparseMatrix/SparseMatrixExample_forAllElements.cpp b/Documentation/Examples/Matrices/SparseMatrix/SparseMatrixExample_forAllElements.cpp
index c603fe32f9..8c9fb368c1 100644
--- a/Documentation/Examples/Matrices/SparseMatrix/SparseMatrixExample_forAllElements.cpp
+++ b/Documentation/Examples/Matrices/SparseMatrix/SparseMatrixExample_forAllElements.cpp
@@ -8,11 +8,9 @@ void forAllElementsExample()
 {
    TNL::Matrices::SparseMatrix< double, Device > matrix( { 1, 2, 3, 4, 5 }, 5 );
 
-   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int& columnIdx, double& value, bool& compute ) {
-      if( rowIdx < localIdx )  // This is important, some matrix formats may allocate more matrix elements
+   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int& columnIdx, double& value ) {
+      if( rowIdx >= localIdx )  // This is important, some matrix formats may allocate more matrix elements
                                // than we requested. These padding elements are processed here as well.
-         compute = false;
-      else
       {
          columnIdx = localIdx;
          value = rowIdx + localIdx + 1;
diff --git a/Documentation/Examples/Matrices/SparseMatrix/SparseMatrixExample_forElements.cpp b/Documentation/Examples/Matrices/SparseMatrix/SparseMatrixExample_forElements.cpp
index 216433b637..2d7bbeba55 100644
--- a/Documentation/Examples/Matrices/SparseMatrix/SparseMatrixExample_forElements.cpp
+++ b/Documentation/Examples/Matrices/SparseMatrix/SparseMatrixExample_forElements.cpp
@@ -8,11 +8,9 @@ void forElementsExample()
 {
    TNL::Matrices::SparseMatrix< double, Device > matrix( { 1, 2, 3, 4, 5 }, 5 );
 
-   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int& columnIdx, double& value, bool& compute ) {
-      if( rowIdx < localIdx )  // This is important, some matrix formats may allocate more matrix elements
+   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int& columnIdx, double& value ) {
+      if( rowIdx >= localIdx )  // This is important, some matrix formats may allocate more matrix elements
                                // than we requested. These padding elements are processed here as well.
-         compute = false;
-      else
       {
          columnIdx = localIdx;
          value = rowIdx + localIdx + 1;
diff --git a/Documentation/Examples/Matrices/SparseMatrix/SparseMatrixViewExample_forAllElements.cpp b/Documentation/Examples/Matrices/SparseMatrix/SparseMatrixViewExample_forAllElements.cpp
index 4000107eb3..79fb7890d5 100644
--- a/Documentation/Examples/Matrices/SparseMatrix/SparseMatrixViewExample_forAllElements.cpp
+++ b/Documentation/Examples/Matrices/SparseMatrix/SparseMatrixViewExample_forAllElements.cpp
@@ -9,11 +9,9 @@ void forAllElementsExample()
    TNL::Matrices::SparseMatrix< double, Device > matrix( { 1, 2, 3, 4, 5 }, 5 );
    auto view = matrix.getView();
 
-   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int& columnIdx, double& value, bool& compute ) {
-      if( rowIdx < localIdx )  // This is important, some matrix formats may allocate more matrix elements
+   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int& columnIdx, double& value ) {
+      if( rowIdx >= localIdx )  // This is important, some matrix formats may allocate more matrix elements
                                // than we requested. These padding elements are processed here as well.
-         compute = false;
-      else
       {
          columnIdx = localIdx;
          value = rowIdx + localIdx + 1;
diff --git a/Documentation/Examples/Matrices/SparseMatrix/SparseMatrixViewExample_forElements.cpp b/Documentation/Examples/Matrices/SparseMatrix/SparseMatrixViewExample_forElements.cpp
index 4ffb2ee834..6e296d3dec 100644
--- a/Documentation/Examples/Matrices/SparseMatrix/SparseMatrixViewExample_forElements.cpp
+++ b/Documentation/Examples/Matrices/SparseMatrix/SparseMatrixViewExample_forElements.cpp
@@ -9,11 +9,9 @@ void forElementsExample()
    TNL::Matrices::SparseMatrix< double, Device > matrix( { 1, 2, 3, 4, 5 }, 5 );
    auto view = matrix.getView();
 
-   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int& columnIdx, double& value, bool& compute ) {
-      if( rowIdx < localIdx )  // This is important, some matrix formats may allocate more matrix elements
+   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int& columnIdx, double& value ) {
+      if( rowIdx >= localIdx )  // This is important, some matrix formats may allocate more matrix elements
                                // than we requested. These padding elements are processed here as well.
-         compute = false;
-      else
       {
          columnIdx = localIdx;
          value = rowIdx + localIdx + 1;
diff --git a/Documentation/Examples/Matrices/TridiagonalMatrix/TridiagonalMatrixExample_forAllElements.cpp b/Documentation/Examples/Matrices/TridiagonalMatrix/TridiagonalMatrixExample_forAllElements.cpp
index c29b439a6b..314cd6e4b8 100644
--- a/Documentation/Examples/Matrices/TridiagonalMatrix/TridiagonalMatrixExample_forAllElements.cpp
+++ b/Documentation/Examples/Matrices/TridiagonalMatrix/TridiagonalMatrixExample_forAllElements.cpp
@@ -20,7 +20,7 @@ void forAllElementsExample()
       5,      // number of matrix rows
       5 );    // number of matrix columns
 
-   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int columnIdx, double& value, bool& compute ) {
+   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int columnIdx, double& value ) {
       /***
        * 'forElements' method iterates only over matrix elements lying on given subdiagonals
        * and so we do not need to check anything. The element value can be expressed
diff --git a/Documentation/Examples/Matrices/TridiagonalMatrix/TridiagonalMatrixExample_forElements.cpp b/Documentation/Examples/Matrices/TridiagonalMatrix/TridiagonalMatrixExample_forElements.cpp
index 243e9468eb..b15a9f5818 100644
--- a/Documentation/Examples/Matrices/TridiagonalMatrix/TridiagonalMatrixExample_forElements.cpp
+++ b/Documentation/Examples/Matrices/TridiagonalMatrix/TridiagonalMatrixExample_forElements.cpp
@@ -20,7 +20,7 @@ void forElementsExample()
       5,      // number of matrix rows
       5 );    // number of matrix columns
 
-   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int columnIdx, double& value, bool& compute ) {
+   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int columnIdx, double& value ) {
       /***
        * 'forElements' method iterates only over matrix elements lying on given subdiagonals
        * and so we do not need to check anything. The element value can be expressed
diff --git a/Documentation/Examples/Matrices/TridiagonalMatrix/TridiagonalMatrixViewExample_forAllElements.cpp b/Documentation/Examples/Matrices/TridiagonalMatrix/TridiagonalMatrixViewExample_forAllElements.cpp
index 0ef4304623..8d90c989ef 100644
--- a/Documentation/Examples/Matrices/TridiagonalMatrix/TridiagonalMatrixViewExample_forAllElements.cpp
+++ b/Documentation/Examples/Matrices/TridiagonalMatrix/TridiagonalMatrixViewExample_forAllElements.cpp
@@ -21,7 +21,7 @@ void forAllElementsExample()
       5 );    // number of matrix columns
    auto view = matrix.getView();
 
-   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int columnIdx, double& value, bool& compute ) {
+   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int columnIdx, double& value ) {
       /***
        * 'forElements' method iterates only over matrix elements lying on given subdiagonals
        * and so we do not need to check anything. The element value can be expressed
diff --git a/Documentation/Examples/Matrices/TridiagonalMatrix/TridiagonalMatrixViewExample_forElements.cpp b/Documentation/Examples/Matrices/TridiagonalMatrix/TridiagonalMatrixViewExample_forElements.cpp
index 3045bc655e..b077c008c7 100644
--- a/Documentation/Examples/Matrices/TridiagonalMatrix/TridiagonalMatrixViewExample_forElements.cpp
+++ b/Documentation/Examples/Matrices/TridiagonalMatrix/TridiagonalMatrixViewExample_forElements.cpp
@@ -21,7 +21,7 @@ void forElementsExample()
       5 );    // number of matrix columns
    auto view = matrix.getView();
 
-   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int columnIdx, double& value, bool& compute ) {
+   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int columnIdx, double& value ) {
       /***
        * 'forElements' method iterates only over matrix elements lying on given subdiagonals
        * and so we do not need to check anything. The element value can be expressed
diff --git a/Documentation/Tutorials/Matrices/DenseMatrixSetup_Benchmark.cpp b/Documentation/Tutorials/Matrices/DenseMatrixSetup_Benchmark.cpp
index 7696e9d0d7..64979b0d41 100644
--- a/Documentation/Tutorials/Matrices/DenseMatrixSetup_Benchmark.cpp
+++ b/Documentation/Tutorials/Matrices/DenseMatrixSetup_Benchmark.cpp
@@ -62,7 +62,7 @@ void forElements( const int matrixSize, Matrix& matrix )
 {
    matrix.setDimensions( matrixSize, matrixSize );
 
-   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int& columnIdx, float& value, bool& compute ) mutable {
+   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int& columnIdx, float& value ) mutable {
       value = rowIdx + columnIdx;
    };
    matrix.forElements( 0, matrixSize, f );
diff --git a/Documentation/Tutorials/Matrices/MultidiagonalMatrixSetup_Benchmark.cpp b/Documentation/Tutorials/Matrices/MultidiagonalMatrixSetup_Benchmark.cpp
index d323105cd0..75186957ea 100644
--- a/Documentation/Tutorials/Matrices/MultidiagonalMatrixSetup_Benchmark.cpp
+++ b/Documentation/Tutorials/Matrices/MultidiagonalMatrixSetup_Benchmark.cpp
@@ -143,7 +143,7 @@ void forElements( const int gridSize, Matrix& matrix )
    const int matrixSize = gridSize * gridSize;
    matrix.setDimensions( matrixSize, matrixSize, getOffsets< typename Matrix::DeviceType >( gridSize ) );
 
-   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int columnIdx, float& value, bool& compute ) mutable {
+   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int columnIdx, float& value ) mutable {
       const int i = rowIdx % gridSize;
       const int j = rowIdx / gridSize;
       if( ( i == 0 || j == 0 || i == gridSize - 1 || j == gridSize - 1 ) && localIdx == 0 )
diff --git a/Documentation/Tutorials/Matrices/SparseMatrixSetup_Benchmark.cpp b/Documentation/Tutorials/Matrices/SparseMatrixSetup_Benchmark.cpp
index 7af7de1e10..31a2a039c0 100644
--- a/Documentation/Tutorials/Matrices/SparseMatrixSetup_Benchmark.cpp
+++ b/Documentation/Tutorials/Matrices/SparseMatrixSetup_Benchmark.cpp
@@ -168,7 +168,7 @@ void forElements( const int gridSize, Matrix& matrix )
    matrix.setDimensions( matrixSize, matrixSize );
    matrix.setRowCapacities( rowCapacities );
 
-   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int& columnIdx, float& value, bool& compute ) mutable {
+   auto f = [=] __cuda_callable__ ( int rowIdx, int localIdx, int& columnIdx, float& value ) mutable {
       const int i = rowIdx % gridSize;
       const int j = rowIdx / gridSize;
       if( ( i == 0 || j == 0 || i == gridSize - 1 || j == gridSize - 1 ) && localIdx == 0 )
diff --git a/src/TNL/Algorithms/Segments/BiEllpackView.h b/src/TNL/Algorithms/Segments/BiEllpackView.h
index 0bb603a3ae..c37ed6d730 100644
--- a/src/TNL/Algorithms/Segments/BiEllpackView.h
+++ b/src/TNL/Algorithms/Segments/BiEllpackView.h
@@ -159,8 +159,7 @@ class BiEllpackView
                 typename Reduction,
                 typename ResultKeeper,
                 typename Real,
-                int BlockDim,
-                typename... Args >
+                int BlockDim >
       __device__
       void reduceSegmentsKernelWithAllParameters( IndexType gridIdx,
                                                      IndexType first,
@@ -168,15 +167,13 @@ class BiEllpackView
                                                      Fetch fetch,
                                                      Reduction reduction,
                                                      ResultKeeper keeper,
-                                                     Real zero,
-                                                     Args... args ) const;
+                                                     Real zero ) const;
 
       template< typename Fetch,
                 typename Reduction,
                 typename ResultKeeper,
                 typename Real_,
-                int BlockDim,
-                typename... Args >
+                int BlockDim >
       __device__
       void reduceSegmentsKernel( IndexType gridIdx,
                                     IndexType first,
@@ -184,8 +181,7 @@ class BiEllpackView
                                     Fetch fetch,
                                     Reduction reduction,
                                     ResultKeeper keeper,
-                                    Real_ zero,
-                                    Args... args ) const;
+                                    Real_ zero ) const;
 
       template< typename View_,
                 typename Index_,
@@ -193,8 +189,7 @@ class BiEllpackView
                 typename Reduction_,
                 typename ResultKeeper_,
                 typename Real_,
-                int BlockDim,
-                typename... Args_ >
+                int BlockDim >
       friend __global__
       void BiEllpackreduceSegmentsKernel( View_ chunkedEllpack,
                                              Index_ gridIdx,
@@ -203,8 +198,7 @@ class BiEllpackView
                                              Fetch_ fetch,
                                              Reduction_ reduction,
                                              ResultKeeper_ keeper,
-                                             Real_ zero,
-                                             Args_... args );
+                                             Real_ zero );
 
       template< typename Index_, typename Fetch_, int BlockDim_, int WarpSize_, bool B_ >
       friend struct details::BiEllpackreduceSegmentsDispatcher;
diff --git a/src/TNL/Algorithms/Segments/BiEllpackView.hpp b/src/TNL/Algorithms/Segments/BiEllpackView.hpp
index 98c5c05c83..b480deac08 100644
--- a/src/TNL/Algorithms/Segments/BiEllpackView.hpp
+++ b/src/TNL/Algorithms/Segments/BiEllpackView.hpp
@@ -275,9 +275,8 @@ forElements( IndexType first, IndexType last, Function&& f ) const
       const IndexType groupsCount = detail::BiEllpack< IndexType, DeviceType, Organization, getWarpSize() >::getActiveGroupsCountDirect( segmentsPermutationView, segmentIdx );
       IndexType groupHeight = getWarpSize();
       //printf( "segmentIdx = %d strip = %d firstGroupInStrip = %d rowStripPerm = %d groupsCount = %d \n", segmentIdx, strip, firstGroupInStrip, rowStripPerm, groupsCount );
-      bool compute( true );
       IndexType localIdx( 0 );
-      for( IndexType groupIdx = firstGroupInStrip; groupIdx < firstGroupInStrip + groupsCount && compute; groupIdx++ )
+      for( IndexType groupIdx = firstGroupInStrip; groupIdx < firstGroupInStrip + groupsCount; groupIdx++ )
       {
          IndexType groupOffset = groupPointersView[ groupIdx ];
          const IndexType groupSize = groupPointersView[ groupIdx + 1 ] - groupOffset;
@@ -289,14 +288,14 @@ forElements( IndexType first, IndexType last, Function&& f ) const
             {
                if( Organization == RowMajorOrder )
                {
-                  f( segmentIdx, localIdx, groupOffset + rowStripPerm * groupWidth + i, compute );
+                  f( segmentIdx, localIdx, groupOffset + rowStripPerm * groupWidth + i );
                }
                else
                {
                   /*printf( "segmentIdx = %d localIdx = %d globalIdx = %d groupIdx = %d groupSize = %d groupWidth = %d\n",
                      segmentIdx, localIdx, groupOffset + rowStripPerm + i * groupHeight,
                      groupIdx, groupSize, groupWidth );*/
-                  f( segmentIdx, localIdx, groupOffset + rowStripPerm + i * groupHeight, compute );
+                  f( segmentIdx, localIdx, groupOffset + rowStripPerm + i * groupHeight );
                }
                localIdx++;
             }
diff --git a/src/TNL/Algorithms/Segments/CSR.h b/src/TNL/Algorithms/Segments/CSR.h
index 4fe24934f2..af05a9f614 100644
--- a/src/TNL/Algorithms/Segments/CSR.h
+++ b/src/TNL/Algorithms/Segments/CSR.h
@@ -293,7 +293,7 @@ class CSR
        */
       OffsetsContainer& getOffsets();
 
-      /***
+      /**
        * \brief Go over all segments and for each segment element call
        * function 'f'. The return type of 'f' is bool.
        * When its true, the for-loop continues. Once 'f' returns false, the for-loop
diff --git a/src/TNL/Algorithms/Segments/CSRView.hpp b/src/TNL/Algorithms/Segments/CSRView.hpp
index 7de1938373..08822ca948 100644
--- a/src/TNL/Algorithms/Segments/CSRView.hpp
+++ b/src/TNL/Algorithms/Segments/CSRView.hpp
@@ -193,9 +193,8 @@ forElements( IndexType begin, IndexType end, Function&& f ) const
       const IndexType begin = offsetsView[ segmentIdx ];
       const IndexType end = offsetsView[ segmentIdx + 1 ];
       IndexType localIdx( 0 );
-      bool compute( true );
-      for( IndexType globalIdx = begin; globalIdx < end && compute; globalIdx++  )
-         f( segmentIdx, localIdx++, globalIdx, compute );
+      for( IndexType globalIdx = begin; globalIdx < end; globalIdx++  )
+         f( segmentIdx, localIdx++, globalIdx );
    };
    Algorithms::ParallelFor< Device >::exec( begin, end, l );
 }
diff --git a/src/TNL/Algorithms/Segments/ChunkedEllpackView.h b/src/TNL/Algorithms/Segments/ChunkedEllpackView.h
index 8f8177dbb4..123bc1cb92 100644
--- a/src/TNL/Algorithms/Segments/ChunkedEllpackView.h
+++ b/src/TNL/Algorithms/Segments/ChunkedEllpackView.h
@@ -158,8 +158,7 @@ class ChunkedEllpackView
       template< typename Fetch,
                 typename Reduction,
                 typename ResultKeeper,
-                typename Real,
-                typename... Args >
+                typename Real >
       __device__
       void reduceSegmentsKernelWithAllParameters( IndexType gridIdx,
                                                      IndexType first,
@@ -167,14 +166,12 @@ class ChunkedEllpackView
                                                      Fetch fetch,
                                                      Reduction reduction,
                                                      ResultKeeper keeper,
-                                                     Real zero,
-                                                     Args... args ) const;
+                                                     Real zero ) const;
 
       template< typename Fetch,
                 typename Reduction,
                 typename ResultKeeper,
-                typename Real,
-                typename... Args >
+                typename Real >
       __device__
       void reduceSegmentsKernel( IndexType gridIdx,
                                     IndexType first,
@@ -182,8 +179,7 @@ class ChunkedEllpackView
                                     Fetch fetch,
                                     Reduction reduction,
                                     ResultKeeper keeper,
-                                    Real zero,
-                                    Args... args ) const;
+                                    Real zero ) const;
 #endif
 
       IndexType size = 0, storageSize = 0, numberOfSlices = 0;
@@ -216,8 +212,7 @@ class ChunkedEllpackView
                 typename Fetch_,
                 typename Reduction_,
                 typename ResultKeeper_,
-                typename Real_,
-                typename... Args_ >
+                typename Real_ >
       friend __global__
       void ChunkedEllpackreduceSegmentsKernel( View_ chunkedEllpack,
                                                   Index_ gridIdx,
@@ -226,8 +221,7 @@ class ChunkedEllpackView
                                                   Fetch_ fetch,
                                                   Reduction_ reduction,
                                                   ResultKeeper_ keeper,
-                                                  Real_ zero,
-                                                  Args_... args );
+                                                  Real_ zero );
 
       template< typename Index_, typename Fetch_, bool B_ >
       friend struct details::ChunkedEllpackreduceSegmentsDispatcher;
diff --git a/src/TNL/Algorithms/Segments/ChunkedEllpackView.hpp b/src/TNL/Algorithms/Segments/ChunkedEllpackView.hpp
index c43b57321d..5f73fd8ab2 100644
--- a/src/TNL/Algorithms/Segments/ChunkedEllpackView.hpp
+++ b/src/TNL/Algorithms/Segments/ChunkedEllpackView.hpp
@@ -323,14 +323,13 @@ forElements( IndexType first, IndexType last, Function&& f ) const
       const IndexType chunkSize = slices[ sliceIdx ].chunkSize;
 
       const IndexType segmentSize = segmentChunksCount * chunkSize;
-      bool compute( true );
       if( Organization == RowMajorOrder )
       {
          IndexType begin = sliceOffset + firstChunkOfSegment * chunkSize;
          IndexType end = begin + segmentSize;
          IndexType localIdx( 0 );
-         for( IndexType j = begin; j < end && compute; j++ )
-            f( segmentIdx, localIdx++, j, compute );
+         for( IndexType j = begin; j < end; j++ )
+            f( segmentIdx, localIdx++, j );
       }
       else
       {
@@ -339,9 +338,9 @@ forElements( IndexType first, IndexType last, Function&& f ) const
          {
             IndexType begin = sliceOffset + firstChunkOfSegment + chunkIdx;
             IndexType end = begin + chunksInSlice * chunkSize;
-            for( IndexType j = begin; j < end && compute; j += chunksInSlice )
+            for( IndexType j = begin; j < end; j += chunksInSlice )
             {
-               f( segmentIdx, localIdx++, j, compute );
+               f( segmentIdx, localIdx++, j );
             }
          }
       }
diff --git a/src/TNL/Algorithms/Segments/EllpackView.hpp b/src/TNL/Algorithms/Segments/EllpackView.hpp
index 120f8f4369..1ec928336b 100644
--- a/src/TNL/Algorithms/Segments/EllpackView.hpp
+++ b/src/TNL/Algorithms/Segments/EllpackView.hpp
@@ -207,9 +207,8 @@ forElements( IndexType first, IndexType last, Function&& f ) const
          const IndexType begin = segmentIdx * segmentSize;
          const IndexType end = begin + segmentSize;
          IndexType localIdx( 0 );
-         bool compute( true );
-         for( IndexType globalIdx = begin; globalIdx < end && compute; globalIdx++  )
-            f( segmentIdx, localIdx++, globalIdx, compute );
+         for( IndexType globalIdx = begin; globalIdx < end; globalIdx++  )
+            f( segmentIdx, localIdx++, globalIdx );
       };
       Algorithms::ParallelFor< Device >::exec( first, last, l );
    }
@@ -221,9 +220,8 @@ forElements( IndexType first, IndexType last, Function&& f ) const
          const IndexType begin = segmentIdx;
          const IndexType end = storageSize;
          IndexType localIdx( 0 );
-         bool compute( true );
-         for( IndexType globalIdx = begin; globalIdx < end && compute; globalIdx += alignedSize )
-            f( segmentIdx, localIdx++, globalIdx, compute );
+         for( IndexType globalIdx = begin; globalIdx < end; globalIdx += alignedSize )
+            f( segmentIdx, localIdx++, globalIdx );
       };
       Algorithms::ParallelFor< Device >::exec( first, last, l );
    }
diff --git a/src/TNL/Algorithms/Segments/SlicedEllpackView.hpp b/src/TNL/Algorithms/Segments/SlicedEllpackView.hpp
index 703da8d3c3..871aa2da05 100644
--- a/src/TNL/Algorithms/Segments/SlicedEllpackView.hpp
+++ b/src/TNL/Algorithms/Segments/SlicedEllpackView.hpp
@@ -242,15 +242,14 @@ forElements( IndexType first, IndexType last, Function&& f ) const
          const IndexType begin = sliceOffsets_view[ sliceIdx ] + segmentInSliceIdx * segmentSize;
          const IndexType end = begin + segmentSize;
          IndexType localIdx( 0 );
-         bool compute( true );
-         for( IndexType globalIdx = begin; globalIdx < end && compute; globalIdx++  )
+         for( IndexType globalIdx = begin; globalIdx < end; globalIdx++  )
          {
             // The following is a workaround of a bug in nvcc 11.2
 #if CUDART_VERSION == 11020
-             f( segmentIdx, localIdx, globalIdx, compute );
+             f( segmentIdx, localIdx, globalIdx );
              localIdx++;
 #else
-             f( segmentIdx, localIdx++, globalIdx, compute );
+             f( segmentIdx, localIdx++, globalIdx );
 #endif
          }
       };
@@ -265,15 +264,14 @@ forElements( IndexType first, IndexType last, Function&& f ) const
          const IndexType begin = sliceOffsets_view[ sliceIdx ] + segmentInSliceIdx;
          const IndexType end = sliceOffsets_view[ sliceIdx + 1 ];
          IndexType localIdx( 0 );
-         bool compute( true );
-         for( IndexType globalIdx = begin; globalIdx < end && compute; globalIdx += SliceSize )
+         for( IndexType globalIdx = begin; globalIdx < end; globalIdx += SliceSize )
          {
             // The following is a workaround of a bug in nvcc 11.2
 #if CUDART_VERSION == 11020
-            f( segmentIdx, localIdx, globalIdx, compute );
+            f( segmentIdx, localIdx, globalIdx );
             localIdx++;
 #else
-            f( segmentIdx, localIdx++, globalIdx, compute );
+            f( segmentIdx, localIdx++, globalIdx );
 #endif
          }
       };
diff --git a/src/TNL/Matrices/DenseMatrix.h b/src/TNL/Matrices/DenseMatrix.h
index a65c12d80b..0036ba240d 100644
--- a/src/TNL/Matrices/DenseMatrix.h
+++ b/src/TNL/Matrices/DenseMatrix.h
@@ -457,10 +457,8 @@ class DenseMatrix : public Matrix< Real, Device, Index, RealAllocator >
        *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
-       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value, bool& compute )`.
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value )`.
        *  The column index repeats twice only for compatibility with sparse matrices.
-       *  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.
@@ -479,10 +477,8 @@ class DenseMatrix : public Matrix< Real, Device, Index, RealAllocator >
        *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
-       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, RealType& value, bool& compute )`.
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, RealType& value )`.
        *  The column index repeats twice only for compatibility with sparse matrices.
-       *  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.
@@ -633,10 +629,8 @@ class DenseMatrix : public Matrix< Real, Device, Index, RealAllocator >
        *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
-       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value, bool& compute )`.
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value )`.
        *  The column index repeats twice only for compatibility with sparse matrices.
-       *  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.
@@ -650,10 +644,8 @@ class DenseMatrix : public Matrix< Real, Device, Index, RealAllocator >
        *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
-       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, RealType& value, bool& compute )`.
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, RealType& value )`.
        *  The column index repeats twice only for compatibility with sparse matrices.
-       *  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.
diff --git a/src/TNL/Matrices/DenseMatrix.hpp b/src/TNL/Matrices/DenseMatrix.hpp
index 16b844cda7..46c85a373b 100644
--- a/src/TNL/Matrices/DenseMatrix.hpp
+++ b/src/TNL/Matrices/DenseMatrix.hpp
@@ -1139,7 +1139,7 @@ operator=( const DenseMatrixView< RHSReal, RHSDevice, RHSIndex, RHSOrganization
    auto this_view = this->view;
    if( std::is_same< DeviceType, RHSDeviceType >::value )
    {
-      auto f = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx, RHSIndexType columnIdx, const RHSRealType& value, bool& compute ) mutable {
+      auto f = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx, RHSIndexType columnIdx, const RHSRealType& value ) mutable {
          this_view( rowIdx, columnIdx ) = value;
       };
       matrix.forAllElements( f );
@@ -1162,7 +1162,7 @@ operator=( const DenseMatrixView< RHSReal, RHSDevice, RHSIndex, RHSOrganization
 
          ////
          // Copy matrix elements into buffer
-         auto f1 = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx, RHSIndexType columnIdx, const RHSRealType& value, bool& compute ) mutable {
+         auto f1 = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx, RHSIndexType columnIdx, const RHSRealType& value ) mutable {
             const IndexType bufferIdx = ( rowIdx - baseRow ) * maxRowLength + columnIdx;
             matrixValuesBuffer_view[ bufferIdx ] = value;
          };
@@ -1214,7 +1214,7 @@ operator=( const RHSMatrix& matrix )
    if( std::is_same< DeviceType, RHSDeviceType >::value )
    {
       const auto segments_view = this->segments.getView();
-      auto f = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx_, RHSIndexType columnIdx, const RHSRealType& value, bool& compute ) mutable {
+      auto f = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx_, RHSIndexType columnIdx, const RHSRealType& value ) mutable {
          if( value != 0.0 && columnIdx != padding_index )
             values_view[ segments_view.getGlobalIndex( rowIdx, columnIdx ) ] = value;
       };
@@ -1244,7 +1244,7 @@ operator=( const RHSMatrix& matrix )
 
          ////
          // Copy matrix elements into buffer
-         auto f1 = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx, RHSIndexType columnIndex, const RHSRealType& value, bool& compute ) mutable {
+         auto f1 = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx, RHSIndexType columnIndex, const RHSRealType& value ) mutable {
             if( columnIndex != padding_index )
             {
                const IndexType bufferIdx = ( rowIdx - baseRow ) * maxRowLength + localIdx;
diff --git a/src/TNL/Matrices/DenseMatrixView.h b/src/TNL/Matrices/DenseMatrixView.h
index ea7f6dbe74..5fd6cda8e3 100644
--- a/src/TNL/Matrices/DenseMatrixView.h
+++ b/src/TNL/Matrices/DenseMatrixView.h
@@ -518,7 +518,7 @@ class DenseMatrixView : public MatrixView< Real, Device, Index >
        *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
-       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx, const RealType& value, bool& compute )`.
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx, const RealType& value )`.
        *  The column index repeats twice only for compatibility with sparse matrices.
        *  If the 'compute' variable is set to false the iteration over the row can
        *  be interrupted.
@@ -540,7 +540,7 @@ class DenseMatrixView : public MatrixView< Real, Device, Index >
        *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
-       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx, RealType& value, bool& compute )`.
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx, RealType& value )`.
        *  The column index repeats twice only for compatibility with sparse matrices.
        *  If the 'compute' variable is set to false the iteration over the row can
        *  be interrupted.
@@ -694,7 +694,7 @@ class DenseMatrixView : public MatrixView< Real, Device, Index >
        *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
-       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value, bool& compute )`.
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value )`.
        *  The column index repeats twice only for compatibility with sparse matrices.
        *  If the 'compute' variable is set to false the iteration over the row can
        *  be interrupted.
@@ -711,7 +711,7 @@ class DenseMatrixView : public MatrixView< Real, Device, Index >
        *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
-       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, RealType& value, bool& compute )`.
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, RealType& value )`.
        *  The column index repeats twice only for compatibility with sparse matrices.
        *  If the 'compute' variable is set to false the iteration over the row can
        *  be interrupted.
diff --git a/src/TNL/Matrices/DenseMatrixView.hpp b/src/TNL/Matrices/DenseMatrixView.hpp
index 68b2de7eed..1ace770568 100644
--- a/src/TNL/Matrices/DenseMatrixView.hpp
+++ b/src/TNL/Matrices/DenseMatrixView.hpp
@@ -355,8 +355,8 @@ DenseMatrixView< Real, Device, Index, Organization >::
 forElements( IndexType begin, IndexType end, Function&& function ) const
 {
    const auto values_view = this->values.getConstView();
-   auto f = [=] __cuda_callable__ ( IndexType rowIdx, IndexType columnIdx, IndexType globalIdx, bool& compute ) mutable {
-      function( rowIdx, columnIdx, columnIdx, values_view[ globalIdx ], compute );
+   auto f = [=] __cuda_callable__ ( IndexType rowIdx, IndexType columnIdx, IndexType globalIdx ) mutable {
+      function( rowIdx, columnIdx, columnIdx, values_view[ globalIdx ] );
    };
    this->segments.forElements( begin, end, f );
 }
@@ -371,8 +371,8 @@ DenseMatrixView< Real, Device, Index, Organization >::
 forElements( IndexType begin, IndexType end, Function&& function )
 {
    auto values_view = this->values.getView();
-   auto f = [=] __cuda_callable__ ( IndexType rowIdx, IndexType columnIdx, IndexType globalIdx, bool& compute ) mutable {
-      function( rowIdx, columnIdx, globalIdx, values_view[ globalIdx ], compute );
+   auto f = [=] __cuda_callable__ ( IndexType rowIdx, IndexType columnIdx, IndexType globalIdx ) mutable {
+      function( rowIdx, columnIdx, globalIdx, values_view[ globalIdx ] );
    };
    this->segments.forElements( begin, end, f );
 }
diff --git a/src/TNL/Matrices/LambdaMatrix.h b/src/TNL/Matrices/LambdaMatrix.h
index 01d3a0b910..cfdd77d2e3 100644
--- a/src/TNL/Matrices/LambdaMatrix.h
+++ b/src/TNL/Matrices/LambdaMatrix.h
@@ -260,10 +260,8 @@ class LambdaMatrix
        *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
-       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value, bool& compute )`.
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value )`.
        *  The column index repeats twice only for compatibility with sparse matrices.
-       *  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.
@@ -348,7 +346,7 @@ class LambdaMatrix
        *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
-       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value, bool& compute )`.
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value )`.
        *  The column index repeats twice only for compatibility with sparse matrices.
        *  If the 'compute' variable is set to false the iteration over the row can
        *  be interrupted.
diff --git a/src/TNL/Matrices/LambdaMatrix.hpp b/src/TNL/Matrices/LambdaMatrix.hpp
index 867016d474..77de0872c5 100644
--- a/src/TNL/Matrices/LambdaMatrix.hpp
+++ b/src/TNL/Matrices/LambdaMatrix.hpp
@@ -317,14 +317,13 @@ forElements( IndexType first, IndexType last, Function& function ) const
    auto matrixElements = this->matrixElementsLambda;
    auto processRow = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
       const IndexType rowLength = rowLengths( rows, columns, rowIdx );
-      bool compute( true );
-      for( IndexType localIdx = 0; localIdx < rowLength && compute; localIdx++ )
+      for( IndexType localIdx = 0; localIdx < rowLength; localIdx++ )
       {
         IndexType elementColumn( 0 );
         RealType elementValue( 0.0 );
         matrixElements( rows, columns, rowIdx, localIdx, elementColumn, elementValue );
         if( elementValue != 0.0 )
-            function( rowIdx, localIdx, elementColumn, elementValue, compute );
+            function( rowIdx, localIdx, elementColumn, elementValue );
       }
    };
    Algorithms::ParallelFor< DeviceType >::exec( first, last, processRow );
diff --git a/src/TNL/Matrices/MultidiagonalMatrix.h b/src/TNL/Matrices/MultidiagonalMatrix.h
index d938a10629..622b79da4f 100644
--- a/src/TNL/Matrices/MultidiagonalMatrix.h
+++ b/src/TNL/Matrices/MultidiagonalMatrix.h
@@ -715,7 +715,7 @@ class MultidiagonalMatrix : public Matrix< Real, Device, Index, RealAllocator >
        * \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 )`,
+       *  `function( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value )`,
        *
        * where
        *
@@ -728,9 +728,6 @@ class MultidiagonalMatrix : public Matrix< Real, Device, Index, RealAllocator >
        *
        * \e value is the matrix element value.
        *
-       * \e compute is a reference to a boolen variable. If it 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.
@@ -749,7 +746,7 @@ class MultidiagonalMatrix : public Matrix< Real, Device, Index, RealAllocator >
        * \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 )`,
+       *  `function( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value )`,
        *
        * where
        *
@@ -762,9 +759,6 @@ class MultidiagonalMatrix : public Matrix< Real, Device, Index, RealAllocator >
        *
        * \e value is a reference to the matrix element value. It can be used even for changing the matrix element value.
        *
-       * \e compute is a reference to a boolen variable. If it 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.
@@ -914,10 +908,8 @@ class MultidiagonalMatrix : public Matrix< Real, Device, Index, RealAllocator >
        *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
-       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value, bool& compute )`.
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value )`.
        *  The column index repeats twice only for compatibility with sparse matrices.
-       *  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.
@@ -931,10 +923,8 @@ class MultidiagonalMatrix : public Matrix< Real, Device, Index, RealAllocator >
        *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
-       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, RealType& value, bool& compute )`.
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, RealType& value )`.
        *  The column index repeats twice only for compatibility with sparse matrices.
-       *  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.
diff --git a/src/TNL/Matrices/MultidiagonalMatrix.hpp b/src/TNL/Matrices/MultidiagonalMatrix.hpp
index b432814c6a..244188831c 100644
--- a/src/TNL/Matrices/MultidiagonalMatrix.hpp
+++ b/src/TNL/Matrices/MultidiagonalMatrix.hpp
@@ -873,7 +873,7 @@ operator=( const MultidiagonalMatrix< Real_, Device_, Index_, Organization_, Rea
       if( std::is_same< Device, Device_ >::value )
       {
          const auto matrix_view = matrix.getView();
-         auto f = [=] __cuda_callable__ ( const IndexType& rowIdx, const IndexType& localIdx, const IndexType& column, Real& value, bool& compute ) mutable {
+         auto f = [=] __cuda_callable__ ( const IndexType& rowIdx, const IndexType& localIdx, const IndexType& column, Real& value ) mutable {
             value = matrix_view.getValues()[ matrix_view.getIndexer().getGlobalIndex( rowIdx, localIdx ) ];
          };
          this->forAllElements( f );
@@ -898,7 +898,7 @@ operator=( const MultidiagonalMatrix< Real_, Device_, Index_, Organization_, Rea
 
             ////
             // Copy matrix elements into buffer
-            auto f1 = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx, RHSIndexType columnIndex, const RHSRealType& value, bool& compute ) mutable {
+            auto f1 = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx, RHSIndexType columnIndex, const RHSRealType& value ) mutable {
                   const IndexType bufferIdx = ( rowIdx - baseRow ) * maxRowLength + localIdx;
                   matrixValuesBuffer_view[ bufferIdx ] = value;
             };
@@ -910,7 +910,7 @@ operator=( const MultidiagonalMatrix< Real_, Device_, Index_, Organization_, Rea
 
             ////
             // Copy matrix elements from the buffer to the matrix
-            auto f2 = [=] __cuda_callable__ ( const IndexType rowIdx, const IndexType localIdx, const IndexType columnIndex, RealType& value, bool& compute  ) mutable {
+            auto f2 = [=] __cuda_callable__ ( const IndexType rowIdx, const IndexType localIdx, const IndexType columnIndex, RealType& value ) mutable {
                const IndexType bufferIdx = ( rowIdx - baseRow ) * maxRowLength + localIdx;
                   value = thisValuesBuffer_view[ bufferIdx ];
             };
diff --git a/src/TNL/Matrices/MultidiagonalMatrixView.h b/src/TNL/Matrices/MultidiagonalMatrixView.h
index 3575602136..fc5799fa68 100644
--- a/src/TNL/Matrices/MultidiagonalMatrixView.h
+++ b/src/TNL/Matrices/MultidiagonalMatrixView.h
@@ -477,7 +477,7 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
        * \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 )`,
+       *  `function( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value )`,
        *
        * where
        *
@@ -490,9 +490,6 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
        *
        * \e value is the matrix element value.
        *
-       * \e compute is a reference to a boolen variable. If it 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.
@@ -511,7 +508,7 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
        * \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 )`,
+       *  `function( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value )`,
        *
        * where
        *
@@ -524,9 +521,6 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
        *
        * \e value is a reference to the matrix element value. It can be used even for changing the matrix element value.
        *
-       * \e compute is a reference to a boolen variable. If it 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.
@@ -676,10 +670,8 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
        *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
-       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value, bool& compute )`.
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value )`.
        *  The column index repeats twice only for compatibility with sparse matrices.
-       *  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.
@@ -693,10 +685,8 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index >
        *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
-       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, RealType& value, bool& compute )`.
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, RealType& value )`.
        *  The column index repeats twice only for compatibility with sparse matrices.
-       *  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.
diff --git a/src/TNL/Matrices/MultidiagonalMatrixView.hpp b/src/TNL/Matrices/MultidiagonalMatrixView.hpp
index 7621f16787..550158e4d4 100644
--- a/src/TNL/Matrices/MultidiagonalMatrixView.hpp
+++ b/src/TNL/Matrices/MultidiagonalMatrixView.hpp
@@ -216,7 +216,7 @@ setValue( const RealType& v )
    // we dont do this->values = v here because it would set even elements 'outside' the matrix
    // method getNumberOfNonzeroElements would not work well then
    const RealType newValue = v;
-   auto f = [=] __cuda_callable__ ( const IndexType& rowIdx, const IndexType& localIdx, const IndexType columnIdx, RealType& value, bool& compute ) mutable {
+   auto f = [=] __cuda_callable__ ( const IndexType& rowIdx, const IndexType& localIdx, const IndexType columnIdx, RealType& value ) mutable {
       value = newValue;
    };
    this->forAllElements( f );
@@ -443,13 +443,12 @@ forElements( IndexType first, IndexType last, Function& function ) const
    const IndexType diagonalsCount = this->diagonalsOffsets.getSize();
    const IndexType columns = this->getColumns();
    const auto indexer = this->indexer;
-   bool compute( true );
    auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
       for( IndexType localIdx = 0; localIdx < diagonalsCount; localIdx++ )
       {
          const IndexType columnIdx = rowIdx + diagonalsOffsets_view[ localIdx ];
          if( columnIdx >= 0 && columnIdx < columns )
-            function( rowIdx, localIdx, columnIdx, values_view[ indexer.getGlobalIndex( rowIdx, localIdx ) ], compute );
+            function( rowIdx, localIdx, columnIdx, values_view[ indexer.getGlobalIndex( rowIdx, localIdx ) ] );
       }
    };
    Algorithms::ParallelFor< DeviceType >::exec( first, last, f );
@@ -469,13 +468,12 @@ forElements( IndexType first, IndexType last, Function& function )
    const IndexType diagonalsCount = this->diagonalsOffsets.getSize();
    const IndexType columns = this->getColumns();
    const auto indexer = this->indexer;
-   bool compute( true );
    auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
-      for( IndexType localIdx = 0; localIdx < diagonalsCount && compute; localIdx++ )
+      for( IndexType localIdx = 0; localIdx < diagonalsCount; localIdx++ )
       {
          const IndexType columnIdx = rowIdx + diagonalsOffsets_view[ localIdx ];
          if( columnIdx >= 0 && columnIdx < columns )
-            function( rowIdx, localIdx, columnIdx, values_view[ indexer.getGlobalIndex( rowIdx, localIdx ) ], compute );
+            function( rowIdx, localIdx, columnIdx, values_view[ indexer.getGlobalIndex( rowIdx, localIdx ) ] );
       }
    };
    Algorithms::ParallelFor< DeviceType >::exec( first, last, f );
diff --git a/src/TNL/Matrices/SparseMatrix.h b/src/TNL/Matrices/SparseMatrix.h
index 237417d66a..ef7fe05803 100644
--- a/src/TNL/Matrices/SparseMatrix.h
+++ b/src/TNL/Matrices/SparseMatrix.h
@@ -724,12 +724,10 @@ class SparseMatrix : public Matrix< Real, Device, Index, RealAllocator >
        * The lambda function `function` should be declared like follows:
        *
        * ```
-       * auto function = [] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value, bool& compute ) { ... };
+       * auto function = [] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value ) { ... };
        * ```
        *
        *  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.
        *
        * \par Example
        * \include Matrices/SparseMatrix/SparseMatrixExample_forElements.cpp
@@ -751,12 +749,10 @@ class SparseMatrix : public Matrix< Real, Device, Index, RealAllocator >
        * The lambda function `function` should be declared like follows:
        *
        * ```
-       * auto function = [] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value, bool& compute ) mutable { ... }
+       * auto function = [] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value ) mutable { ... }
        * ```
        *
        *  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.
        *
        * \par Example
        * \include Matrices/SparseMatrix/SparseMatrixExample_forElements.cpp
@@ -903,10 +899,8 @@ class SparseMatrix : public Matrix< Real, Device, Index, RealAllocator >
        *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
-       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value, bool& compute )`.
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value )`.
        *  The column index repeats twice only for compatibility with sparse matrices.
-       *  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.
@@ -920,10 +914,8 @@ class SparseMatrix : public Matrix< Real, Device, Index, RealAllocator >
        *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
-       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, RealType& value, bool& compute )`.
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, RealType& value )`.
        *  The column index repeats twice only for compatibility with sparse matrices.
-       *  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.
diff --git a/src/TNL/Matrices/SparseMatrix.hpp b/src/TNL/Matrices/SparseMatrix.hpp
index 08e55190f4..40c568b4df 100644
--- a/src/TNL/Matrices/SparseMatrix.hpp
+++ b/src/TNL/Matrices/SparseMatrix.hpp
@@ -890,7 +890,7 @@ operator=( const DenseMatrix< Real_, Device_, Index_, Organization, RealAllocato
    if( std::is_same< DeviceType, RHSDeviceType >::value )
    {
       const auto segments_view = this->segments.getView();
-      auto f = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx, RHSIndexType columnIdx, const RHSRealType& value, bool& compute ) mutable {
+      auto f = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx, RHSIndexType columnIdx, const RHSRealType& value ) mutable {
          if( value != 0.0 )
          {
             IndexType thisGlobalIdx = segments_view.getGlobalIndex( rowIdx, rowLocalIndexes_view[ rowIdx ]++ );
@@ -921,7 +921,7 @@ operator=( const DenseMatrix< Real_, Device_, Index_, Organization, RealAllocato
 
          ////
          // Copy matrix elements into buffer
-         auto f1 = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx, RHSIndexType columnIndex, const RHSRealType& value, bool& compute ) mutable {
+         auto f1 = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx, RHSIndexType columnIndex, const RHSRealType& value ) mutable {
             const IndexType bufferIdx = ( rowIdx - baseRow ) * maxRowLength + localIdx;
             matrixValuesBuffer_view[ bufferIdx ] = value;
          };
@@ -935,7 +935,7 @@ operator=( const DenseMatrix< Real_, Device_, Index_, Organization, RealAllocato
          // Copy matrix elements from the buffer to the matrix and ignoring
          // zero matrix elements.
          const IndexType matrix_columns = this->getColumns();
-         auto f2 = [=] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, IndexType& columnIndex, RealType& value, bool& compute  ) mutable {
+         auto f2 = [=] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, IndexType& columnIndex, RealType& value ) mutable {
             RealType inValue( 0.0 );
             IndexType bufferIdx, column( rowLocalIndexes_view[ rowIdx ] );
             while( inValue == 0.0 && column < matrix_columns )
@@ -1001,7 +1001,7 @@ operator=( const RHSMatrix& matrix )
    if( std::is_same< DeviceType, RHSDeviceType >::value )
    {
       const auto segments_view = this->segments.getView();
-      auto f = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx_, RHSIndexType columnIndex, const RHSRealType& value, bool& compute ) mutable {
+      auto f = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx_, RHSIndexType columnIndex, const RHSRealType& value ) mutable {
          IndexType localIdx( rowLocalIndexes_view[ rowIdx ] );
          if( value != 0.0 && columnIndex != paddingIndex )
          {
@@ -1043,7 +1043,7 @@ operator=( const RHSMatrix& matrix )
 
          ////
          // Copy matrix elements into buffer
-         auto f1 = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx, RHSIndexType columnIndex, const RHSRealType& value, bool& compute ) mutable {
+         auto f1 = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx, RHSIndexType columnIndex, const RHSRealType& value ) mutable {
             if( columnIndex != paddingIndex )
             {
                TNL_ASSERT_LT( rowIdx - baseRow, bufferRowsCount, "" );
@@ -1066,7 +1066,7 @@ operator=( const RHSMatrix& matrix )
          // zero matrix elements
          //const IndexType matrix_columns = this->getColumns();
          const auto thisRowLengths_view = thisRowLengths.getConstView();
-         auto f2 = [=] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, IndexType& columnIndex, RealType& value, bool& compute ) mutable {
+         auto f2 = [=] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, IndexType& columnIndex, RealType& value ) mutable {
             RealType inValue( 0.0 );
             size_t bufferIdx;
             IndexType bufferLocalIdx( rowLocalIndexes_view[ rowIdx ] );
diff --git a/src/TNL/Matrices/SparseMatrixView.h b/src/TNL/Matrices/SparseMatrixView.h
index 40a89b628a..10310d8023 100644
--- a/src/TNL/Matrices/SparseMatrixView.h
+++ b/src/TNL/Matrices/SparseMatrixView.h
@@ -508,10 +508,8 @@ class SparseMatrixView : public MatrixView< Real, Device, Index >
        *
        * \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 )`.
+       *  `function( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value )`.
        *  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.
@@ -530,10 +528,8 @@ class SparseMatrixView : public MatrixView< Real, Device, Index >
        *
        * \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 )`.
+       *  `function( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value )`.
        *  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.
@@ -685,10 +681,8 @@ class SparseMatrixView : public MatrixView< Real, Device, Index >
        *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
-       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value, bool& compute )`.
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value )`.
        *  The column index repeats twice only for compatibility with sparse matrices.
-       *  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.
@@ -702,10 +696,8 @@ class SparseMatrixView : public MatrixView< Real, Device, Index >
        *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
-       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, RealType& value, bool& compute )`.
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, RealType& value )`.
        *  The column index repeats twice only for compatibility with sparse matrices.
-       *  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.
diff --git a/src/TNL/Matrices/SparseMatrixView.hpp b/src/TNL/Matrices/SparseMatrixView.hpp
index 5df969867f..72377f8478 100644
--- a/src/TNL/Matrices/SparseMatrixView.hpp
+++ b/src/TNL/Matrices/SparseMatrixView.hpp
@@ -594,11 +594,11 @@ forElements( IndexType begin, IndexType end, Function& function ) const
    const auto columns_view = this->columnIndexes.getConstView();
    const auto values_view = this->values.getConstView();
    //const IndexType paddingIndex_ = this->getPaddingIndex();
-   auto f = [=] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, IndexType globalIdx, bool& compute ) mutable -> bool {
+   auto f = [=] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, IndexType globalIdx ) mutable -> bool {
       if( isBinary() )
-         function( rowIdx, localIdx, columns_view[ globalIdx ], 1, compute );
+         function( rowIdx, localIdx, columns_view[ globalIdx ], 1 );
       else
-         function( rowIdx, localIdx, columns_view[ globalIdx ], values_view[ globalIdx ], compute );
+         function( rowIdx, localIdx, columns_view[ globalIdx ], values_view[ globalIdx ] );
       return true;
    };
    this->segments.forElements( begin, end, f );
@@ -618,14 +618,14 @@ forElements( IndexType begin, IndexType end, Function& function )
    auto columns_view = this->columnIndexes.getView();
    auto values_view = this->values.getView();
    const IndexType paddingIndex_ = this->getPaddingIndex();
-   auto f = [=] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, IndexType globalIdx, bool& compute ) mutable {
+   auto f = [=] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, IndexType globalIdx ) mutable {
       if( isBinary() )
       {
          RealType one( columns_view[ globalIdx ] != paddingIndex_ );
-         function( rowIdx, localIdx, columns_view[ globalIdx ], one, compute );
+         function( rowIdx, localIdx, columns_view[ globalIdx ], one );
       }
       else
-         function( rowIdx, localIdx, columns_view[ globalIdx ], values_view[ globalIdx ], compute );
+         function( rowIdx, localIdx, columns_view[ globalIdx ], values_view[ globalIdx ] );
    };
    this->segments.forElements( begin, end, f );
 }
diff --git a/src/TNL/Matrices/TridiagonalMatrix.h b/src/TNL/Matrices/TridiagonalMatrix.h
index b74e0dcb9d..4789a079f0 100644
--- a/src/TNL/Matrices/TridiagonalMatrix.h
+++ b/src/TNL/Matrices/TridiagonalMatrix.h
@@ -610,10 +610,8 @@ class TridiagonalMatrix : public Matrix< Real, Device, Index, RealAllocator >
        *
        * \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 )`.
+       *  `function( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value )`.
        *  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.
@@ -632,10 +630,8 @@ class TridiagonalMatrix : public Matrix< Real, Device, Index, RealAllocator >
        *
        * \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 )`.
+       *  `function( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value )`.
        *  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.
@@ -654,10 +650,8 @@ class TridiagonalMatrix : public Matrix< Real, Device, Index, RealAllocator >
        *
        * \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 )`.
+       *  `function( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value )`.
        *  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.
@@ -676,10 +670,8 @@ class TridiagonalMatrix : public Matrix< Real, Device, Index, RealAllocator >
        *
        * \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 )`.
+       *  `function( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value )`.
        *  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.
@@ -798,10 +790,8 @@ class TridiagonalMatrix : public Matrix< Real, Device, Index, RealAllocator >
        *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
-       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value, bool& compute )`.
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value )`.
        *  The column index repeats twice only for compatibility with sparse matrices.
-       *  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.
@@ -815,10 +805,8 @@ class TridiagonalMatrix : public Matrix< Real, Device, Index, RealAllocator >
        *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
-       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, RealType& value, bool& compute )`.
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, RealType& value )`.
        *  The column index repeats twice only for compatibility with sparse matrices.
-       *  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.
diff --git a/src/TNL/Matrices/TridiagonalMatrix.hpp b/src/TNL/Matrices/TridiagonalMatrix.hpp
index a76ae3ef92..93ebcd8b39 100644
--- a/src/TNL/Matrices/TridiagonalMatrix.hpp
+++ b/src/TNL/Matrices/TridiagonalMatrix.hpp
@@ -717,7 +717,7 @@ operator=( const TridiagonalMatrix< Real_, Device_, Index_, Organization_, RealA
       if( std::is_same< Device, Device_ >::value )
       {
          const auto matrix_view = matrix.getView();
-         auto f = [=] __cuda_callable__ ( const IndexType& rowIdx, const IndexType& localIdx, const IndexType& column, Real& value, bool& compute ) mutable {
+         auto f = [=] __cuda_callable__ ( const IndexType& rowIdx, const IndexType& localIdx, const IndexType& column, Real& value ) mutable {
             value = matrix_view.getValues()[ matrix_view.getIndexer().getGlobalIndex( rowIdx, localIdx ) ];
          };
          this->forAllElements( f );
@@ -727,7 +727,7 @@ operator=( const TridiagonalMatrix< Real_, Device_, Index_, Organization_, RealA
          TridiagonalMatrix< Real, Device, Index, Organization_ > auxMatrix;
          auxMatrix = matrix;
          const auto matrix_view = auxMatrix.getView();
-         auto f = [=] __cuda_callable__ ( const IndexType& rowIdx, const IndexType& localIdx, const IndexType& column, Real& value, bool& compute ) mutable {
+         auto f = [=] __cuda_callable__ ( const IndexType& rowIdx, const IndexType& localIdx, const IndexType& column, Real& value ) mutable {
             value = matrix_view.getValues()[ matrix_view.getIndexer().getGlobalIndex( rowIdx, localIdx ) ];
          };
          this->forAllElements( f );
diff --git a/src/TNL/Matrices/TridiagonalMatrixView.h b/src/TNL/Matrices/TridiagonalMatrixView.h
index e05a8b0597..6b68f41978 100644
--- a/src/TNL/Matrices/TridiagonalMatrixView.h
+++ b/src/TNL/Matrices/TridiagonalMatrixView.h
@@ -463,10 +463,8 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
        *
        * \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 )`.
+       *  `function( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value )`.
        *  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.
@@ -485,10 +483,8 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
        *
        * \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 )`.
+       *  `function( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value )`.
        *  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.
@@ -639,7 +635,7 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
        *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
-       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value, bool& compute )`.
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value )`.
        *  The column index repeats twice only for compatibility with sparse matrices.
        *  If the 'compute' variable is set to false the iteration over the row can
        *  be interrupted.
@@ -656,7 +652,7 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
        *
        * \tparam Function is type of lambda function that will operate on matrix elements.
        *    It is should have form like
-       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, RealType& value, bool& compute )`.
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, RealType& value )`.
        *  The column index repeats twice only for compatibility with sparse matrices.
        *  If the 'compute' variable is set to false the iteration over the row can
        *  be interrupted.
diff --git a/src/TNL/Matrices/TridiagonalMatrixView.hpp b/src/TNL/Matrices/TridiagonalMatrixView.hpp
index 0ebad8cb64..c670733810 100644
--- a/src/TNL/Matrices/TridiagonalMatrixView.hpp
+++ b/src/TNL/Matrices/TridiagonalMatrixView.hpp
@@ -393,26 +393,25 @@ forElements( IndexType first, IndexType last, Function& function ) const
 {
    const auto values_view = this->values.getConstView();
    const auto indexer = this->indexer;
-   bool compute( true );
    auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
       if( rowIdx == 0 )
       {
-         function( 0, 1, 0, values_view[ indexer.getGlobalIndex( 0, 1 ) ], compute );
-         function( 0, 2, 1, values_view[ indexer.getGlobalIndex( 0, 2 ) ], compute );
+         function( 0, 1, 0, values_view[ indexer.getGlobalIndex( 0, 1 ) ] );
+         function( 0, 2, 1, values_view[ indexer.getGlobalIndex( 0, 2 ) ] );
       }
       else if( rowIdx + 1 < indexer.getColumns() )
       {
-         function( rowIdx, 0, rowIdx - 1, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ], compute );
-         function( rowIdx, 1, rowIdx,     values_view[ indexer.getGlobalIndex( rowIdx, 1 ) ], compute );
-         function( rowIdx, 2, rowIdx + 1, values_view[ indexer.getGlobalIndex( rowIdx, 2 ) ], compute );
+         function( rowIdx, 0, rowIdx - 1, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ] );
+         function( rowIdx, 1, rowIdx,     values_view[ indexer.getGlobalIndex( rowIdx, 1 ) ] );
+         function( rowIdx, 2, rowIdx + 1, values_view[ indexer.getGlobalIndex( rowIdx, 2 ) ] );
       }
       else if( rowIdx < indexer.getColumns() )
       {
-         function( rowIdx, 0, rowIdx - 1, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ], compute );
-         function( rowIdx, 1, rowIdx,     values_view[ indexer.getGlobalIndex( rowIdx, 1 ) ], compute );
+         function( rowIdx, 0, rowIdx - 1, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ] );
+         function( rowIdx, 1, rowIdx,     values_view[ indexer.getGlobalIndex( rowIdx, 1 ) ] );
       }
       else
-         function( rowIdx, 0, rowIdx, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ], compute );
+         function( rowIdx, 0, rowIdx, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ] );
    };
    Algorithms::ParallelFor< DeviceType >::exec( first, last, f );
 }
@@ -428,26 +427,25 @@ forElements( IndexType first, IndexType last, Function& function )
 {
    auto values_view = this->values.getView();
    const auto indexer = this->indexer;
-   bool compute( true );
    auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
       if( rowIdx == 0 )
       {
-         function( 0, 1, 0, values_view[ indexer.getGlobalIndex( 0, 1 ) ], compute );
-         function( 0, 2, 1, values_view[ indexer.getGlobalIndex( 0, 2 ) ], compute );
+         function( 0, 1, 0, values_view[ indexer.getGlobalIndex( 0, 1 ) ] );
+         function( 0, 2, 1, values_view[ indexer.getGlobalIndex( 0, 2 ) ] );
       }
       else if( rowIdx + 1 < indexer.getColumns() )
       {
-         function( rowIdx, 0, rowIdx - 1, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ], compute );
-         function( rowIdx, 1, rowIdx,     values_view[ indexer.getGlobalIndex( rowIdx, 1 ) ], compute );
-         function( rowIdx, 2, rowIdx + 1, values_view[ indexer.getGlobalIndex( rowIdx, 2 ) ], compute );
+         function( rowIdx, 0, rowIdx - 1, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ] );
+         function( rowIdx, 1, rowIdx,     values_view[ indexer.getGlobalIndex( rowIdx, 1 ) ] );
+         function( rowIdx, 2, rowIdx + 1, values_view[ indexer.getGlobalIndex( rowIdx, 2 ) ] );
       }
       else if( rowIdx < indexer.getColumns() )
       {
-         function( rowIdx, 0, rowIdx - 1, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ], compute );
-         function( rowIdx, 1, rowIdx,     values_view[ indexer.getGlobalIndex( rowIdx, 1 ) ], compute );
+         function( rowIdx, 0, rowIdx - 1, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ] );
+         function( rowIdx, 1, rowIdx,     values_view[ indexer.getGlobalIndex( rowIdx, 1 ) ] );
       }
       else
-         function( rowIdx, 0, rowIdx, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ], compute );
+         function( rowIdx, 0, rowIdx, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ] );
    };
    Algorithms::ParallelFor< DeviceType >::exec( first, last, f );
 }
@@ -663,13 +661,13 @@ addMatrix( const TridiagonalMatrixView< Real_, Device_, Index_, Organization_ >&
       const auto matrix_view = matrix;
       const auto matrixMult = matrixMultiplicator;
       const auto thisMult = thisMatrixMultiplicator;
-      auto add0 = [=] __cuda_callable__ ( const IndexType& rowIdx, const IndexType& localIdx, const IndexType& column, Real& value, bool& compute ) mutable {
+      auto add0 = [=] __cuda_callable__ ( const IndexType& rowIdx, const IndexType& localIdx, const IndexType& column, Real& value ) mutable {
          value = matrixMult * matrix.getValues()[ matrix.getIndexer().getGlobalIndex( rowIdx, localIdx ) ];
       };
-      auto add1 = [=] __cuda_callable__ ( const IndexType& rowIdx, const IndexType& localIdx, const IndexType& column, Real& value, bool& compute ) mutable {
+      auto add1 = [=] __cuda_callable__ ( const IndexType& rowIdx, const IndexType& localIdx, const IndexType& column, Real& value ) mutable {
          value += matrixMult * matrix.getValues()[ matrix.getIndexer().getGlobalIndex( rowIdx, localIdx ) ];
       };
-      auto addGen = [=] __cuda_callable__ ( const IndexType& rowIdx, const IndexType& localIdx, const IndexType& column, Real& value, bool& compute ) mutable {
+      auto addGen = [=] __cuda_callable__ ( const IndexType& rowIdx, const IndexType& localIdx, const IndexType& column, Real& value ) mutable {
          value = thisMult * value + matrixMult * matrix.getValues()[ matrix.getIndexer().getGlobalIndex( rowIdx, localIdx ) ];
       };
       if( thisMult == 0.0 )
diff --git a/src/UnitTests/Matrices/DenseMatrixTest.h b/src/UnitTests/Matrices/DenseMatrixTest.h
index ceb7ae3588..9a5fee6e98 100644
--- a/src/UnitTests/Matrices/DenseMatrixTest.h
+++ b/src/UnitTests/Matrices/DenseMatrixTest.h
@@ -791,7 +791,7 @@ void test_ForElements()
    const IndexType rows = 8;
 
    Matrix m( rows, cols  );
-   m.forAllElements( [] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, const IndexType& columnIdx, RealType& value, bool compute ) mutable {
+   m.forAllElements( [] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, const IndexType& columnIdx, RealType& value ) mutable {
       value = rowIdx + 1.0;
    } );
 
diff --git a/src/UnitTests/Matrices/SparseMatrixTest.hpp b/src/UnitTests/Matrices/SparseMatrixTest.hpp
index cca22d8574..1716b0ab89 100644
--- a/src/UnitTests/Matrices/SparseMatrixTest.hpp
+++ b/src/UnitTests/Matrices/SparseMatrixTest.hpp
@@ -1050,7 +1050,7 @@ void test_ForElements()
    const IndexType rows = 8;
 
    Matrix m( { 3, 3, 3, 3, 3, 3, 3, 3, 3 }, cols  );
-   m.forAllElements( [] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, IndexType& columnIdx, RealType& value, bool compute ) mutable {
+   m.forAllElements( [] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, IndexType& columnIdx, RealType& value ) mutable {
       value = rowIdx + 1.0;
       columnIdx = localIdx;
    } );
diff --git a/src/UnitTests/Matrices/SparseMatrixVectorProductTest.hpp b/src/UnitTests/Matrices/SparseMatrixVectorProductTest.hpp
index c3f8a39dbb..d39593e232 100644
--- a/src/UnitTests/Matrices/SparseMatrixVectorProductTest.hpp
+++ b/src/UnitTests/Matrices/SparseMatrixVectorProductTest.hpp
@@ -365,7 +365,7 @@ void test_VectorProduct_largeMatrix()
       TNL::Containers::Vector< IndexType, DeviceType, IndexType > rowCapacities( size );
       rowCapacities.forAllElements( [] __cuda_callable__ ( IndexType i, IndexType& value ) { value = 1; } );
       m1.setRowCapacities( rowCapacities );
-      auto f1 = [=] __cuda_callable__ ( IndexType row, IndexType localIdx, IndexType& column, RealType& value, bool& compute ) {
+      auto f1 = [=] __cuda_callable__ ( IndexType row, IndexType localIdx, IndexType& column, RealType& value ) {
          if( localIdx == 0  )
          {
             value = row + 1;
@@ -389,7 +389,7 @@ void test_VectorProduct_largeMatrix()
       rowCapacities.setSize( rows );
       rowCapacities.forAllElements( [=] __cuda_callable__ ( IndexType i, IndexType& value ) { value = i + 1; } );
       m2.setRowCapacities( rowCapacities );
-      auto f2 = [=] __cuda_callable__ ( IndexType row, IndexType localIdx, IndexType& column, RealType& value, bool& compute ) {
+      auto f2 = [=] __cuda_callable__ ( IndexType row, IndexType localIdx, IndexType& column, RealType& value ) {
          if( localIdx <= row )
          {
             value = localIdx + 1;
@@ -436,7 +436,7 @@ void test_VectorProduct_longRowsMatrix()
          TNL::Containers::Vector< IndexType, DeviceType, IndexType > rowsCapacities( rows );
          rowsCapacities = columns;
          m3.setRowCapacities( rowsCapacities );
-         auto f = [] __cuda_callable__ ( IndexType row, IndexType localIdx, IndexType& column, RealType& value, bool& compute ) {
+         auto f = [] __cuda_callable__ ( IndexType row, IndexType localIdx, IndexType& column, RealType& value ) {
             column = localIdx;
             value = localIdx + row;
          };
-- 
GitLab