diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 6ebb5cd0c513bd0b45ae740ed922fd023922cce2..9d49ddc6f3254e3ab166e93183b80d9da0f5aff5 100644 --- a/.gitlab-ci.yml +++ b/.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] diff --git a/dune/CMakeLists.txt b/dune/CMakeLists.txt index 5b54675ecb3902e440b371a4fbc25513584a713a..e9087bc607541d4e902724e630b38d257dc0f075 100644 --- a/dune/CMakeLists.txt +++ b/dune/CMakeLists.txt @@ -1,2 +1,5 @@ add_subdirectory(istl) -add_subdirectory(python) +# if Python bindings are enabled, include necessary sub directories. +if( DUNE_ENABLE_PYTHONBINDINGS ) + add_subdirectory("python") +endif() diff --git a/dune/python/CMakeLists.txt b/dune/python/CMakeLists.txt index e93d9cba6bc17643f35959eaac7fb9c7afe603a3..864ed72ca6489c85e76ba56aafafe8ea6cd7b70e 100644 --- a/dune/python/CMakeLists.txt +++ b/dune/python/CMakeLists.txt @@ -1 +1,2 @@ add_subdirectory(istl) +add_subdirectory(test) diff --git a/dune/python/istl/bcrsmatrix.hh b/dune/python/istl/bcrsmatrix.hh index 528d5cf78c1e250e6ec002d4e9b782925b87905d..fcf180c18172efc178422eb1d93a87c33f02d889 100644 --- a/dune/python/istl/bcrsmatrix.hh +++ b/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 diff --git a/dune/python/istl/bvector.hh b/dune/python/istl/bvector.hh index cdbb82b52986a9b7bd6f3a39149d860bddcd3b4e..5f7006de9b1971e0707bbd803de0231aca2b4530 100644 --- a/dune/python/istl/bvector.hh +++ b/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 ); diff --git a/dune/python/istl/matrixindexset.hh b/dune/python/istl/matrixindexset.hh index 71e9f495a99c6eaec165d9ecc1cedf640dfa45c0..d4c81e1be801172fd79a1176fd7a2d11bb23ffba 100644 --- a/dune/python/istl/matrixindexset.hh +++ b/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 diff --git a/dune/python/istl/operators.hh b/dune/python/istl/operators.hh index cb19a10c04c8eba24a17157c23910fb3a4936135..c043da4aa0ae3a775fe50c22816fa0dcaf9bd887 100644 --- a/dune/python/istl/operators.hh +++ b/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> diff --git a/dune/python/istl/preconditioners.hh b/dune/python/istl/preconditioners.hh index b4620e1f369cd6419f647c70abdaecbec21ae1ee..87de453c2ceee4719786f6ddb2ceb01680b03716 100644 --- a/dune/python/istl/preconditioners.hh +++ b/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 diff --git a/dune/python/istl/solvers.hh b/dune/python/istl/solvers.hh index e37ae86d908acc439e39f5da4a3af38b92471585..b4706ec3164155396c640860bda1a252beebcf1c 100644 --- a/dune/python/istl/solvers.hh +++ b/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 diff --git a/dune/python/test/CMakeLists.txt b/dune/python/test/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..65bbdcb51579c0a848add00bf4715759a24ecd48 --- /dev/null +++ b/dune/python/test/CMakeLists.txt @@ -0,0 +1,4 @@ +dune_python_add_test(NAME pybcrsmatrix + COMMAND ${PYTHON_EXECUTABLE} bcrsmatrix.py + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} + LABELS quick) diff --git a/dune/python/test/bcrsmatrix.mm b/dune/python/test/bcrsmatrix.mm new file mode 100644 index 0000000000000000000000000000000000000000..04709ab6b4433323f06d1f80a4e9c1f13349735e --- /dev/null +++ b/dune/python/test/bcrsmatrix.mm @@ -0,0 +1,8 @@ +%%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 diff --git a/dune/python/test/bcrsmatrix.py b/dune/python/test/bcrsmatrix.py new file mode 100644 index 0000000000000000000000000000000000000000..8027f3ff8bc6cb03fd957139e59f74fb46a109a6 --- /dev/null +++ b/dune/python/test/bcrsmatrix.py @@ -0,0 +1,116 @@ +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) diff --git a/dune/python/test/bcrsmatrix.txt b/dune/python/test/bcrsmatrix.txt new file mode 100644 index 0000000000000000000000000000000000000000..0c88019fba71695dc1c146b0c74c109481f8be95 --- /dev/null +++ b/dune/python/test/bcrsmatrix.txt @@ -0,0 +1,5 @@ +1 1 1 +2 2 1 +3 3 1 +4 4 1 +5 5 1