added a test for the BCRSMatrix and some cleanup

.gitlab-ci.yml

@@ -14,3 +14,22 @@ variables:
 # number of cores. We solve the issue by disallowing OpenMP to allocate more
 # than one thread.
   OMP_NUM_THREADS: 1
+
+debian-11-gcc-9-17-python:
+  image: duneci/debian:11
+  script: duneci-standard-test
+  variables:
+    DUNECI_TOOLCHAIN:   gcc-9-17
+    # so we need some variables to build the dune-py module during execution of the first python test:
+    # we need to find the dune mdoule
+    DUNE_CONTROL_PATH:  /duneci/modules:$CI_PROJECT_DIR
+    # the position for the dune-py module
+    DUNE_PY_DIR:        /duneci/modules/dune-py
+    # during dune-py build this variable is used - do know a way to access
+    # the CMAKE_FLAGS used to build the modules...
+    DUNE_CMAKE_FLAGS:   "CC=gcc-9 CXX=g++-9 -DCMAKE_CXX_FLAGS='-std=c++17 -O2 -g -Wall -fdiagnostics-color=always' -DDUNE_ENABLE_PYTHONBINDINGS=ON -DDUNE_MAX_TEST_CORES=4 -DBUILD_SHARED_LIBS=TRUE -DDUNE_PYTHON_INSTALL_LOCATION=none -DCMAKE_POSITION_INDEPENDENT_CODE=TRUE -DCMAKE_DISABLE_FIND_PACKAGE_LATEX=TRUE -DCMAKE_DISABLE_FIND_PACKAGE_Alberta=TRUE -DCMAKE_DISABLE_FIND_PACKAGE_Vc=TRUE -DCMAKE_DISABLE_DOCUMENTATION=TRUE"
+    # cmake flags we use for all dune moudle - important is that build shared libs is set (need some better way of doing this)
+    DUNECI_CMAKE_FLAGS: $DUNE_CMAKE_FLAGS
+    # finally set the python path to all modules
+    PYTHONPATH:         /duneci/modules/dune-common/build-cmake/python:/duneci/modules/dune-geometry/build-cmake/python:$CI_PROJECT_DIR/build-cmake/python
+  tags: [duneci]

dune/CMakeLists.txt

 add_subdirectory(istl)
-add_subdirectory(python)
+# if Python bindings are enabled, include necessary sub directories.
+if( DUNE_ENABLE_PYTHONBINDINGS )
+  add_subdirectory("python")
+endif()

dune/python/CMakeLists.txt

 add_subdirectory(istl)
+add_subdirectory(test)

dune/python/istl/bcrsmatrix.hh

@@ -6,17 +6,14 @@
 #include <string>
 #include <tuple>
 
-#include <dune/common/fmatrix.hh>
-#include <dune/common/fvector.hh>
+#include <dune/python/common/fvecmatregistry.hh>
 
-#if HAVE_DUNE_ISTL
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 #include <dune/istl/io.hh>
 #include <dune/istl/matrixmarket.hh>
 #include <dune/istl/matrixindexset.hh>
 #include <dune/istl/operators.hh>
-#endif // #if HAVE_DUNE_ISTL
 
 #include <dune/python/pybind11/pybind11.h>
 #include <dune/python/pybind11/stl.h>
@@ -74,7 +71,6 @@ namespace Dune
     // registerBCRSMatrix
     // ------------------
 
-#if HAVE_DUNE_ISTL
     template <class BCRSMatrix, class... options>
     void registerBCRSMatrix(pybind11::handle scope,
                             pybind11::class_<BCRSMatrix, options...> cls)
@@ -82,13 +78,12 @@ namespace Dune
       using pybind11::operator""_a;
       typedef typename BCRSMatrix::block_type block_type;
       typedef typename BCRSMatrix::field_type field_type;
-      //typedef typename FieldTraits< field_type >::real_type real_type;
       typedef typename BCRSMatrix::size_type Size;
       typedef typename BCRSMatrix::row_type row_type;
 
       typedef typename BCRSMatrix::BuildMode BuildMode;
-      //need one of thses for istl buildmode
-      //typedef typename istl::BuildMode istlBuildMode;
+
+      registerFieldVecMat<block_type>::apply();
 
       pybind11::class_< row_type > clsRow( scope, "BCRSMatrixRow" );
       registerBlockVector( clsRow );
@@ -180,7 +175,7 @@ namespace Dune
           if( row >= self.N() )
             throw pybind11::index_error( "No such row: " + std::to_string( row ) );
           return self[ row ];
-        }, pybind11::return_value_policy::reference, pybind11::keep_alive< 0, 1 >() );
+        }, pybind11::return_value_policy::reference); // , pybind11::keep_alive< 0, 1 >() );
       cls.def( "__getitem__", [] ( BCRSMatrix &self, std::tuple< Size, Size > index ) -> block_type & {
           const Size row = std::get< 0 >( index );
           if( row >= self.N() )
@@ -340,11 +335,10 @@ namespace Dune
       auto cls = Dune::Python::insertClass< BCRSMatrix >( scope, clsName, Dune::Python::GenerateTypeName(matrixTypename), Dune::Python::IncludeFiles{"dune/istl/bcrsmatrix.hh","dune/python/istl/bcrsmatrix.hh"});
       if(cls.second)
       {
-      registerBCRSMatrix( scope, cls.first );
+        registerBCRSMatrix( scope, cls.first );
       }
       return cls.first;
     }
-#endif // #if HAVE_DUNE_ISTL
 
   } // namespace Python
 

dune/python/istl/bvector.hh

@@ -12,17 +12,14 @@
 
 //this added otherwise insert class wasn't possible on line ~190
 #include <dune/python/common/typeregistry.hh>
-
-
+#include <dune/python/common/fvecmatregistry.hh>
 #include <dune/python/common/string.hh>
 #include <dune/python/common/vector.hh>
 #include <dune/python/istl/iterator.hh>
 #include <dune/python/pybind11/operators.h>
 #include <dune/python/pybind11/pybind11.h>
 
-#if HAVE_DUNE_ISTL
 #include <dune/istl/bvector.hh>
-#endif // #if HAVE_DUNE_ISTL
 
 namespace Dune
 {
@@ -44,7 +41,6 @@ namespace Dune
       }
 
 
-#if HAVE_DUNE_ISTL
       template< class B, class A >
       inline static void copy ( const char *ptr, const ssize_t *shape, const ssize_t *strides, Dune::BlockVector< B, A > &v )
       {
@@ -52,7 +48,6 @@ namespace Dune
         for( ssize_t i = 0; i < *shape; ++i )
           copy( ptr + i*(*strides), shape+1, strides+1, v[ i ] );
       }
-#endif // #if HAVE_DUNE_ISTL
 
 
       template< class BlockVector >
@@ -112,6 +107,8 @@ namespace Dune
       typedef typename BlockVector::block_type block_type;
       typedef typename BlockVector::size_type size_type;
 
+      registerFieldVecMat<block_type>::apply();
+
       using pybind11::operator""_a;
 
       cls.def( "assign", [] ( BlockVector &self, const BlockVector &x ) { self = x; }, "x"_a );

dune/python/istl/matrixindexset.hh

@@ -9,9 +9,7 @@
 #include <dune/common/fmatrix.hh>
 #include <dune/common/fvector.hh>
 
-#if HAVE_DUNE_ISTL
 #include <dune/istl/matrixindexset.hh>
-#endif // #if HAVE_DUNE_ISTL
 
 #include <dune/python/pybind11/pybind11.h>
 #include <dune/python/pybind11/stl.h>
@@ -25,7 +23,6 @@ namespace Dune
     // registermatrixindexset
     // ------------------
 
-#if HAVE_DUNE_ISTL
     template <class MatrixIndexSet, class... options>
     void registerMatrixIndexSet(pybind11::handle scope,
                             pybind11::class_<MatrixIndexSet, options...> cls)
@@ -48,7 +45,6 @@ namespace Dune
       registerMatrixIndexSet( scope, cls );
       return cls;
     }
-#endif // #if HAVE_DUNE_ISTL
 
   } // namespace Python
 

dune/python/istl/operators.hh

@@ -3,9 +3,7 @@
 
 #include <dune/common/typeutilities.hh>
 
-#if HAVE_DUNE_ISTL
 #include <dune/istl/operators.hh>
-#endif // #if HAVE_DUNE_ISTL
 
 #include <dune/python/pybind11/pybind11.h>
 

dune/python/istl/preconditioners.hh

@@ -4,11 +4,9 @@
 #include <dune/common/typeutilities.hh>
 #include <dune/common/version.hh>
 
-#if HAVE_DUNE_ISTL
 #include <dune/istl/operators.hh>
 #include <dune/istl/preconditioner.hh>
 #include <dune/istl/preconditioners.hh>
-#endif // #if HAVE_DUNE_ISTL
 
 #include <dune/python/pybind11/pybind11.h>
 
@@ -21,7 +19,6 @@ namespace Dune
     // registerPreconditioner
     // ----------------------
 
-#if HAVE_DUNE_ISTL
     template< class Preconditioner, class... options >
     inline void registerPreconditioner ( pybind11::class_< Preconditioner, options... > cls )
     {
@@ -47,14 +44,12 @@ namespace Dune
 
       cls.def_property_readonly( "category", [] ( const Preconditioner &self ) { return self.category(); } );
     }
-#endif // #if HAVE_DUNE_ISTL
 
 
 
     // registerPreconditioners
     // -----------------------
 
-#if HAVE_DUNE_ISTL
     template< class X, class Y, class... options >
     inline void registerPreconditioners ( pybind11::module module, pybind11::class_< LinearOperator< X, Y >, options...  > cls )
     {
@@ -84,9 +79,7 @@ namespace Dune
               This is a sequential preconditioner.
         )doc" );
     }
-#endif // #if HAVE_DUNE_ISTL
 
-#if HAVE_DUNE_ISTL
     template< class M, class X, class Y, class... options >
     inline void registerMatrixPreconditioners ( pybind11::module module, pybind11::class_< LinearOperator< X, Y >, options... > cls )
     {
@@ -208,7 +201,6 @@ namespace Dune
               given matrix, however only L and D are actually stored.
         )doc" );
     }
-#endif // #if HAVE_DUNE_ISTL
 
   } // namespace Python
 

dune/python/istl/solvers.hh

@@ -4,11 +4,9 @@
 #include <dune/common/typeutilities.hh>
 #include <dune/common/version.hh>
 
-#if HAVE_DUNE_ISTL
 #include <dune/istl/solver.hh>
 #include <dune/istl/solvers.hh>
 #include <dune/istl/preconditioners.hh>
-#endif // #if HAVE_DUNE_ISTL
 
 #include <dune/python/istl/preconditioners.hh>
 
@@ -84,7 +82,6 @@ namespace Dune
       // registerEndomorphismSolvers
       // ---------------------------
 
-#if HAVE_DUNE_ISTL
       template< class X, class Y, class... options >
       inline std::enable_if_t< std::is_same< X, Y >::value >
       registerEndomorphismSolvers ( pybind11::module module, pybind11::class_< LinearOperator< X, Y >, options... >, PriorityTag< 1 > )
@@ -204,7 +201,6 @@ namespace Dune
       {
         registerEndomorphismSolvers( module, cls, PriorityTag< 42 >() );
       }
-#endif // #if HAVE_DUNE_ISTL
 
     } // namespace detail
 
@@ -213,7 +209,6 @@ namespace Dune
     // registerSolvers
     // ---------------
 
-#if HAVE_DUNE_ISTL
     template< class X, class Y, class... options >
     inline void registerSolvers ( pybind11::module module, pybind11::class_< LinearOperator< X, Y >, options... > cls )
     {
@@ -251,7 +246,6 @@ namespace Dune
                 This can lead to a large memory consumption.
           )doc" );
     }
-#endif // #if HAVE_DUNE_ISTL
 
   } // namespace Python
 

dune/python/test/CMakeLists.txt

+dune_python_add_test(NAME pybcrsmatrix
+                     COMMAND ${PYTHON_EXECUTABLE} bcrsmatrix.py
+                     WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
+                     LABELS quick)

dune/python/test/bcrsmatrix.mm

+%%MatrixMarket matrix coordinate real general
+% ISTL_STRUCT blocked 1 1
+5 5 5
+1 1 1.000000e+00
+2 2 1.000000e+00
+3 3 1.000000e+00
+4 4 1.000000e+00
+5 5 1.000000e+00

dune/python/test/bcrsmatrix.py

+import dune.common
+from dune.istl import blockVector, BlockVector, bcrsMatrix, BCRSMatrix, SeqJacobi, CGSolver, BuildMode
+
+def identity(n, buildMode=BuildMode.random):
+    if buildMode is BuildMode.random:
+        mat = bcrsMatrix((n, n), n, BuildMode.random)
+        for i in range(n):
+            mat.setRowSize(i, 1)
+        mat.endRowSizes()
+        for i in range(n):
+            mat.setIndices(i, [i])
+        mat.endIndices()
+        for i in range(n):
+            mat[i, i] = 1
+        return mat
+    elif buildMode is BuildMode.implicit:
+        mat = bcrsMatrix((n, n), n, 0, BuildMode.implicit)
+        for i in range(n):
+            mat[i, i] = 1
+        mat.compress()
+        return mat
+    else:
+        raise Exception("buildMode " + str(buildMode) + " not supported by identity")
+
+def isIdentity(mat):
+    if mat.rows != mat.cols:
+        return False
+    for i in range(mat.rows):
+        if not (i, i) in mat:
+            return False
+        for j in range(mat.cols):
+            if (i, j) in mat and mat[i, j] != (1 if i == j else 0):
+                return False
+    return True
+
+mat = bcrsMatrix((5, 5), 5, BuildMode.random)
+if mat.shape != (5, 5) or mat.nonZeroes != 5 or mat.buildMode != BuildMode.random:
+    raise Exception("Matrix does not return constructor arguments correctly")
+
+mat = bcrsMatrix((5, 5), 5, BuildMode.implicit, blockType=(2,3))
+if mat.shape != (5, 5) or mat.nonZeroes != 5 or mat.buildMode != BuildMode.implicit \
+   or mat[0,0].rows != 2 or mat[0,0].cols != 3:
+    raise Exception("Matrix does not return constructor arguments correctly")
+
+Matrix = BCRSMatrix((2,3))
+mat = Matrix((5, 5), 5, Matrix.BuildMode.implicit)
+if mat.shape != (5, 5) or mat.nonZeroes != 5 or mat.buildMode != BuildMode.implicit \
+   or mat[0,0].rows != 2 or mat[0,0].cols != 3:
+    raise Exception("Matrix does not return constructor arguments correctly")
+
+# the following not working yet since the matrix market exported fails
+# BMatrix = BCRSMatrix(Matrix)
+# mat = BMatrix((5,5), 5, BuildMode.implicit)
+
+# store identity matrix
+mat = identity(5)
+if not isIdentity(mat):
+    raise Exception("Identity matrix not setup correctly")
+mat.store("bcrsmatrix.mm", "matrixmarket")
+
+for i, row in mat.enumerate:
+    for j, col in row.enumerate:
+        if i != j:
+            raise Exception("Wrong sparsity pattern")
+        if col != 1:
+            raise Exception("Diagonal entry is not 1")
+
+mat = identity(5, BuildMode.implicit)
+if not isIdentity(mat):
+    raise Exception("Identity matrix not setup correctly")
+
+# manipulate diagonal to 2
+for i in range(mat.rows):
+    mat[i, i] *= 2
+if isIdentity(mat):
+    raise Exception("Matrix not manipulated")
+
+# reload identity matrix
+mat = bcrsMatrix()
+mat.load("bcrsmatrix.mm", "matrixmarket")
+if not isIdentity(mat):
+    raise Exception("Matrix not loaded correctly")
+
+# store in matlab format
+mat.store("bcrsmatrix.txt", "matlab")
+
+# solve trivial linear system
+x = blockVector(5)
+for i in range(5):
+    x[i] = i+1 # float(i+1) # dune.common.FieldVector([i+1])
+
+y1, y2 = blockVector(5), blockVector(5)
+mat.mv(x, y1)
+mat.asLinearOperator().apply(x, y2)
+if (y1 - y2).two_norm > 1e-12:
+    raise Exception("mat.mv != mat.asLinearOperator().apply")
+
+S = CGSolver(mat.asLinearOperator(), SeqJacobi(mat), 1e-10)
+mat = None
+z = blockVector(5)
+_, _, converged, _, _ = S(z, y1)
+if not converged:
+    raise Exception("CGSolver has not converged")
+if (z - x).two_norm > 1e-8:
+    raise Exception("CGSolver unable to solve identity")
+
+s = "(" + ", ".join(str(v) for v in x) + ")"
+
+str_x = "("
+for i in range(0,5):
+    str_x = str_x + "(" +  "{0:.6f}".format(x[i][0]) + "), "
+str_x = str_x[:-2]
+str_x = str_x +")"
+
+if str_x != s:
+    raise Exception(str(x) + " = str(x) != " + s)

dune/python/test/bcrsmatrix.txt

+1 1 1
+2 2 1
+3 3 1
+4 4 1
+5 5 1