diff --git a/Documentation/Tutorials/Matrices/BinarySparseMatrixExample.cpp b/Documentation/Tutorials/Matrices/BinarySparseMatrixExample.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ef17f7044216740416b62cfe581f9513c8807b77
--- /dev/null
+++ b/Documentation/Tutorials/Matrices/BinarySparseMatrixExample.cpp
@@ -0,0 +1,41 @@
+#include <iostream>
+#include <TNL/Matrices/SparseMatrix.h>
+#include <TNL/Algorithms/Segments/CSR.h>
+#include <TNL/Devices/Host.h>
+
+
+template< typename Device >
+void binarySparseMatrixExample()
+{
+   TNL::Matrices::SparseMatrix< bool, Device, int > binaryMatrix (
+      5, // number of matrix rows
+      5, // number of matrix columns
+      {  // matrix elements definition
+         {  0,  0, 1.0 }, {  0,  1, 2.0 }, {  0,  2, 3.0 }, {  0,  3, 4.0 }, {  0,  4, 5.0 },
+         {  1,  0, 2.0 }, {  1,  1,  1.0 },
+         {  2,  0, 3.0 }, {  2,  2,  1.0 },
+         {  3,  0, 4.0 }, {  3,  3,  1.0 },
+         {  4,  0, 5.0 }, {  4,  4,  1.0 } } );
+
+   std::cout << "Binary sparse matrix: " << std::endl << binaryMatrix << std::endl;
+
+   TNL::Containers::Vector< double, Device > inVector( 5, 1.1 ), outVector( 5, 0.0 );
+   binaryMatrix.vectorProduct( inVector, outVector );
+   std::cout << "Product with vector " << inVector << " is " << outVector << std::endl << std::endl;
+
+   TNL::Matrices::SparseMatrix< bool, Device, int, TNL::Matrices::GeneralMatrix, TNL::Algorithms::Segments::CSRDefault, double > binaryMatrix2;
+   binaryMatrix2 = binaryMatrix;
+   binaryMatrix2.vectorProduct( inVector, outVector );
+   std::cout << "Product with vector in double precision " << inVector << " is " << outVector << std::endl << std::endl;
+}
+
+int main( int argc, char* argv[] )
+{
+   std::cout << "Creating matrix on CPU ... " << std::endl;
+   binarySparseMatrixExample< TNL::Devices::Host >();
+
+#ifdef HAVE_CUDA
+   std::cout << "Creating matrix on CUDA GPU ... " << std::endl;
+   binarySparseMatrixExample< TNL::Devices::Cuda >();
+#endif
+}
diff --git a/Documentation/Tutorials/Matrices/BinarySparseMatrixExample.cu b/Documentation/Tutorials/Matrices/BinarySparseMatrixExample.cu
new file mode 120000
index 0000000000000000000000000000000000000000..4311752ce7d46340befbb1b9732f26facd1df771
--- /dev/null
+++ b/Documentation/Tutorials/Matrices/BinarySparseMatrixExample.cu
@@ -0,0 +1 @@
+BinarySparseMatrixExample.cpp
\ No newline at end of file
diff --git a/Documentation/Tutorials/Matrices/CMakeLists.txt b/Documentation/Tutorials/Matrices/CMakeLists.txt
index 176ade630a157af49dec26a08e429c000fb296c2..9a48b5afa25ea2053cc7f75422af21e9294a9f61 100644
--- a/Documentation/Tutorials/Matrices/CMakeLists.txt
+++ b/Documentation/Tutorials/Matrices/CMakeLists.txt
@@ -98,6 +98,13 @@ IF( BUILD_CUDA )
    ADD_CUSTOM_COMMAND( COMMAND SymmetricSparseMatrixExample >
                         ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/SymmetricSparseMatrixExample.out
                         OUTPUT SymmetricSparseMatrixExample.out )
+
+   CUDA_ADD_EXECUTABLE( BinarySparseMatrixExample BinarySparseMatrixExample.cu )
+   ADD_CUSTOM_COMMAND( COMMAND BinarySparseMatrixExample >
+                        ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/BinarySparseMatrixExample.out
+                        OUTPUT BinarySparseMatrixExample.out )
+
+
    ####
    # THe following examples/benchmarks run for very long time
    CUDA_ADD_EXECUTABLE( DenseMatrixSetup_Benchmark_cuda DenseMatrixSetup_Benchmark.cu )
@@ -133,6 +140,7 @@ ADD_CUSTOM_TARGET( TutorialsMatricesCuda ALL DEPENDS
    SparseMatrixExample_rowsReduction_vectorProduct.out
    SparseMatrixViewExample_setElement.out
    SymmetricSparseMatrixExample.out
+   BinarySparseMatrixExample.out
  )
 ELSE()
 ADD_CUSTOM_TARGET( TutorialsMatrices ALL DEPENDS
diff --git a/Documentation/Tutorials/Matrices/tutorial_Matrices.md b/Documentation/Tutorials/Matrices/tutorial_Matrices.md
index 5f3c01d737d0025cb132be56c73fb909ceaf02a8..df9d38258c9cdd3d3248f2d188ff0ba87d2a1342 100644
--- a/Documentation/Tutorials/Matrices/tutorial_Matrices.md
+++ b/Documentation/Tutorials/Matrices/tutorial_Matrices.md
@@ -100,7 +100,7 @@ There is no change in the dense matrix part of the table. The numbers grow propo
 | Real   | Index  | Dense matrix | Multidiagonal matrix |  Sparse matrix | Fill ratio |
 |:------:|:------:|:------------:|:--------------------:|:--------------:|:----------:|
 | float  | 32-bit |          4 B |                  4 B |            8 B |     << 50% |
-| float  | 32-bit |          4 B |                  4 B |           12 B |     << 30% |
+| float  | 64-bit |          4 B |                  4 B |           12 B |     << 30% |
 | double | 32-bit |          8 B |                  8 B |           12 B |     << 60% |
 | double | 64-bit |          8 B |                  8 B |           16 B |     << 50% |
 
@@ -632,7 +632,7 @@ would not make sense. If we pass through this test, the matrix element lies in t
 
 #### Symmetric sparse matrices
 
-For sparse [symmetric matrices](https://en.wikipedia.org/wiki/Symmetric_matrix), TNL offers a format storing only a half of the matrix elements. More precisely, ony the matrix diagonal and the elements bellow are stored in the memory. The matrix elements above the diagonal are deduced from those bellow. If such a symmetric format is used on GPU, atomic operations must be used in some matrix operations. For this reason, symmetric matrices are allowed only for when the matrix elements values are expressed with `float` and `double` type. An advantage of the symmetric formats is lower memory consumption. Since less data need to be transferred from the memory, better performance might be observed. In some cases, however, the use of atomic operations on GPU may cause performance drop. Mostly we can see approximately the same performance compared to general formats but we can profit from lower memory requirements which is appreciated especially on GPU. The following example shows how to create symmetric sparse matrix.
+For sparse [symmetric matrices](https://en.wikipedia.org/wiki/Symmetric_matrix), TNL offers a format storing only a half of the matrix elements. More precisely, ony the matrix diagonal and the elements bellow are stored in the memory. The matrix elements above the diagonal are deduced from those bellow. If such a symmetric format is used on GPU, atomic operations must be used in some matrix operations. For this reason, symmetric matrices can be combined only with matrix elements values expressed in `float` or `double` type. An advantage of the symmetric formats is lower memory consumption. Since less data need to be transferred from the memory, better performance might be observed. In some cases, however, the use of atomic operations on GPU may cause performance drop. Mostly we can see approximately the same performance compared to general formats but we can profit from lower memory requirements which is appreciated especially on GPU. The following example shows how to create symmetric sparse matrix.
 
 \includelineno SymmetricSparseMatrixExample.cpp
 
@@ -656,6 +656,24 @@ The elements depicted in grey color are not stored in the memory. The main diffe
 
 **Warning: Assignment of symmetric sparse matrix to general sparse matrix does not give correct result, currently. Only the diagonal and the lower part of the matrix is assigned.**
 
+#### Binary sparse matrices
+
+If the matrix element value type (i.e. `Real` type) is set to `bool` the matrix elements can be only `1` or `0`. So in the sparse matrix formats, where we do not store the zero matrix elements, explicitly stored elements can have only one possible value which is `1`.  Therefore we do not need to store the values, only the positions of the nonzero elements. The array `values`, which usualy stores the matrix elements values, can be completely omitted and we can reduce the memory requirements. The following table shows how much we can reduce the memory consumption when using binary matrix instead of common sparse matrix using `float` or `double` types:
+
+| Real   | Index  | Common sparse matrix | Binary sparse matrix | Ratio      |
+|:------:|:------:|:--------------------:|:--------------------:|:----------:|
+| float  | 32-bit |         4 + 4 =  8 B |                  4 B |        50% |
+| float  | 64-bit |         4 + 8 = 12 B |                  8 B |        75% |
+| double | 32-bit |         8 + 4 = 12 B |                  4 B |        33% |
+| double | 64-bit |         8 + 8 = 16 B |                  8 B |        50% |
+
+The following example demonstrates the use of binary matrix:
+
+\includelineno BinarySparseMatrixExample.cpp
+
+All we need to do is set the `Real` type to `bool` as we can see on the line 9. We can see that even though we set different values to different matrix elements (lines 14-18) at the end all of them are turned into ones (printing of the matrix on the line 20). There is an issue, however, which is demonstrated on the product of the matrix with a vector. Nonbinary matrices compute all operations using the `Real` type. If it is set to `bool` operations like [SpMV](https://en.wikipedia.org/wiki/Sparse_matrix-vector_multiplication) would not get correct solution. Therefore sparse matrices use another type called `ComputeReal` which is the 6th template parameter of \ref TNL::Matrices::SparseMatrix. By default it is set to `Index` type but it can be changed by the user. On the lines 26-29 we show how to change this type to `double` and what is the effect of it (correct result of matrix-vector multiplication). The result looks as follows:
+
+\include BinarySparseMatrixExample.out
 
 ### Tridiagonal matrices <a name="tridiagonal_matrices_setup"></a>
 
diff --git a/src/TNL/Matrices/SparseMatrix.h b/src/TNL/Matrices/SparseMatrix.h
index 044c4698b25b40cc1ef8e411f85692279928cbf5..2100b05a3d5d0a06045014a28df555f24c13db00 100644
--- a/src/TNL/Matrices/SparseMatrix.h
+++ b/src/TNL/Matrices/SparseMatrix.h
@@ -139,14 +139,14 @@ class SparseMatrix : public Matrix< Real, Device, Index, RealAllocator >
        *
        * See \ref SparseMatrixView.
        */
-      using ViewType = SparseMatrixView< Real, Device, Index, MatrixType, SegmentsViewTemplate >;
+      using ViewType = SparseMatrixView< Real, Device, Index, MatrixType, SegmentsViewTemplate, ComputeRealType >;
 
       /**
        * \brief Matrix view type for constant instances.
        *
        * See \ref SparseMatrixView.
        */
-      using ConstViewType = SparseMatrixView< std::add_const_t< Real >, Device, Index, MatrixType, SegmentsViewTemplate >;
+      using ConstViewType = SparseMatrixView< std::add_const_t< Real >, Device, Index, MatrixType, SegmentsViewTemplate, ComputeRealType >;
 
       /**
        * \brief Type for accessing matrix rows.
diff --git a/src/TNL/Matrices/SparseMatrixRowView.h b/src/TNL/Matrices/SparseMatrixRowView.h
index e54c8bf89b95b163809478856c7f2666b7807783..84da4e064c44c3198c7665cb2334d7a0ae1d0efa 100644
--- a/src/TNL/Matrices/SparseMatrixRowView.h
+++ b/src/TNL/Matrices/SparseMatrixRowView.h
@@ -19,19 +19,19 @@ namespace Matrices {
 
 /**
  * \brief RowView is a simple structure for accessing rows of sparse matrix.
- * 
+ *
  * \tparam SegmentView is a segment view of segments representing the matrix format.
  * \tparam ValuesView is a vector view storing the matrix elements values.
  * \tparam ColumnsIndexesView is a vector view storing the column indexes of the matrix element.
  * \tparam isBinary tells if the the parent matrix is a binary matrix.
- * 
+ *
  * See \ref SparseMatrix and \ref SparseMatrixView.
- * 
+ *
  * \par Example
  * \include Matrices/SparseMatrix/SparseMatrixExample_getRow.cpp
  * \par Output
  * \include SparseMatrixExample_getRow.out
- * 
+ *
  * \par Example
  * \include Matrices/SparseMatrix/SparseMatrixViewExample_getRow.cpp
  * \par Output
@@ -87,13 +87,13 @@ class SparseMatrixRowView
 
       /**
        * \brief Tells whether the parent matrix is a binary matrix.
-       * @return 
+       * @return `true` if the matrix is binary.
        */
       static constexpr bool isBinary() { return isBinary_; };
 
       /**
        * \brief Constructor with \e segmentView, \e values and \e columnIndexes.
-       * 
+       *
        * \param segmentView instance of SegmentViewType representing matrix row.
        * \param values is a container view for storing the matrix elements values.
        * \param columnIndexes is a container view for storing the column indexes of the matrix elements.
@@ -105,7 +105,7 @@ class SparseMatrixRowView
 
       /**
        * \brief Returns size of the matrix row, i.e. number of matrix elements in this row.
-       * 
+       *
        * \return Size of the matrix row.
        */
       __cuda_callable__
@@ -113,9 +113,9 @@ class SparseMatrixRowView
 
       /**
        * \brief Returns constants reference to a column index of an element with given rank in the row.
-       * 
+       *
        * \param localIdx is the rank of the non-zero element in given row.
-       * 
+       *
        * \return constant reference to the matrix element column index.
        */
       __cuda_callable__
@@ -123,9 +123,9 @@ class SparseMatrixRowView
 
       /**
        * \brief Returns non-constants reference to a column index of an element with given rank in the row.
-       * 
+       *
        * \param localIdx is the rank of the non-zero element in given row.
-       * 
+       *
        * \return non-constant reference to the matrix element column index.
        */
       __cuda_callable__
@@ -133,9 +133,9 @@ class SparseMatrixRowView
 
       /**
        * \brief Returns constants reference to value of an element with given rank in the row.
-       * 
+       *
        * \param localIdx is the rank of the non-zero element in given row.
-       * 
+       *
        * \return constant reference to the matrix element value.
        */
       __cuda_callable__
@@ -143,9 +143,9 @@ class SparseMatrixRowView
 
       /**
        * \brief Returns non-constants reference to value of an element with given rank in the row.
-       * 
+       *
        * \param localIdx is the rank of the non-zero element in given row.
-       * 
+       *
        * \return non-constant reference to the matrix element value.
        */
       __cuda_callable__
@@ -153,7 +153,7 @@ class SparseMatrixRowView
 
       /**
        * \brief Sets a value of matrix element with given rank in the matrix row.
-       * 
+       *
        * \param localIdx is the rank of the matrix element in the row.
        * \param value is the new value of the matrix element.
        */
@@ -163,7 +163,7 @@ class SparseMatrixRowView
 
       /**
        * \brief Sets a column index of matrix element with given rank in the matrix row.
-       * 
+       *
        * \param localIdx is the rank of the matrix element in the row.
        * \param columnIndex is the new column index of the matrix element.
        */
@@ -173,7 +173,7 @@ class SparseMatrixRowView
 
       /**
        * \brief Sets both a value and a column index of matrix element with given rank in the matrix row.
-       * 
+       *
        * \param localIdx is the rank of the matrix element in the row.
        * \param columnIndex is the new column index of the matrix element.
        * \param value is the new value of the matrix element.
@@ -185,9 +185,9 @@ class SparseMatrixRowView
 
       /**
        * \brief Comparison of two matrix rows.
-       * 
+       *
        * The other matrix row can be from any other matrix.
-       * 
+       *
        * \param other is another matrix row.
        * \return \e true if both rows are the same, \e false otherwise.
        */
@@ -209,7 +209,7 @@ class SparseMatrixRowView
 
 /**
  * \brief Insertion operator for a sparse matrix row.
- * 
+ *
  * \param str is an output stream.
  * \param row is an input sparse matrix row.
  * \return  reference to the output stream.