diff --git a/src/FiniteDifference/Derivative.hpp b/src/FiniteDifference/Derivative.hpp index b356276fddb2333f154ed75ff53ae70b14584e4a..c3795f63844910e56bbef389ca16b82d8a3a9316 100644 --- a/src/FiniteDifference/Derivative.hpp +++ b/src/FiniteDifference/Derivative.hpp @@ -18,7 +18,7 @@ #include "Grid/comb.hpp" #include "FiniteDifference/util/common.hpp" #include "util/util_num.hpp" - +#include <unordered_map> /*! \brief Derivative second order on h (spacing) * diff --git a/src/FiniteDifference/FDScheme.hpp b/src/FiniteDifference/FDScheme.hpp index 9fe06bdece05b48a996db1dff5abaedd4874d3c5..06dbcf55bae762ec9895d917634450a99d703a48 100644 --- a/src/FiniteDifference/FDScheme.hpp +++ b/src/FiniteDifference/FDScheme.hpp @@ -15,6 +15,9 @@ #include "data_type/scalar.hpp" #include "NN/CellList/CellDecomposer.hpp" #include "Grid/staggered_dist_grid_util.hpp" +#include "Grid/grid_dist_id.hpp" +#include "Vector/Vector_util.hpp" +#include "Grid/staggered_dist_grid.hpp" /*! \brief Finite Differences * diff --git a/src/Makefile.am b/src/Makefile.am index 867a536aed8e748d3d1f5dbe006592289bf87936..2f0fdc1beae5e7e654cc536e92c4193a8d471d41 100755 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,11 +1,11 @@ -LINKLIBS = $(PETSC_LIB) $(SUITESPARSE_LIBS) $(LAPACK_LIBS) $(BLAS_LIBS) $(METIS_LIB) $(PARMETIS_LIB) $(DEFAULT_LIB) $(PTHREAD_LIBS) $(OPT_LIBS) $(BOOST_LDFLAGS) $(BOOST_PROGRAM_OPTIONS_LIB) $(BOOST_IOSTREAMS_LIB) $(HDF5_LDFLAGS) $(HDF5_LIBS) $(LIBQUADMATH) $(PETSC_LIB) +LINKLIBS = $(PETSC_LIB) $(SUITESPARSE_LIBS) $(LAPACK_LIBS) $(BLAS_LIBS) $(METIS_LIB) $(PARMETIS_LIB) $(DEFAULT_LIB) $(PTHREAD_LIBS) $(OPT_LIBS) $(BOOST_LDFLAGS) $(BOOST_PROGRAM_OPTIONS_LIB) $(BOOST_IOSTREAMS_LIB) $(HDF5_LDFLAGS) $(HDF5_LIBS) $(LIBQUADMATH) noinst_PROGRAMS = numerics numerics_SOURCES = main.cpp ../../openfpm_vcluster/src/VCluster.cpp ../../openfpm_devices/src/memory/HeapMemory.cpp ../../openfpm_devices/src/memory/PtrMemory.cpp ../../openfpm_devices/src/Memleak_check.cpp -numerics_CXXFLAGS = $(AM_CXXFLAGS) $(HDF5_CPPFLAGS) $(INCLUDES_PATH) $(BOOST_CPPFLAGS) $(SUITESPARSE_INCLUDE) $(METIS_INCLUDE) $(PARMETIS_INCLUDE) $(EIGEN_INCLUDE) $(PETSC_INCLUDE) -Wno-deprecated-declarations -Wno-unused-local-typedefs +numerics_CXXFLAGS = $(AM_CXXFLAGS) $(HDF5_CPPFLAGS) $(INCLUDES_PATH) $(BOOST_CPPFLAGS) $(SUITESPARSE_INCLUDE) $(METIS_INCLUDE) $(PARMETIS_INCLUDE) $(EIGEN_INCLUDE) $(PETSC_INCLUDE) -Wno-deprecated-declarations -Wno-unused-local-typedefs numerics_CFLAGS = $(CUDA_CFLAGS) -numerics_LDADD = $(LINKLIBS) +numerics_LDADD = $(LINKLIBS) -lmetis -lparmetis nobase_include_HEADERS = PSE/Kernels.hpp PSE/Kernels_test_util.hpp .cu.o : diff --git a/src/Matrix/SparseMatrix.hpp b/src/Matrix/SparseMatrix.hpp index ca60651c4be5117bdd85bf34aa45892707ca4a44..6b067bcce17817fb82e8a50e6e37edbe0f8e18a8 100644 --- a/src/Matrix/SparseMatrix.hpp +++ b/src/Matrix/SparseMatrix.hpp @@ -8,7 +8,12 @@ #ifndef OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_HPP_ #define OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_HPP_ +#ifdef HAVE_EIGEN #include <Eigen/Sparse> +#define DEFAULT_MATRIX = Eigen::SparseMatrix<T,0,id_t> +#else +#define DEFAULT_MATRIX = void +#endif /*! \brief It store the non zero elements of the matrix * @@ -30,6 +35,17 @@ template<typename T, unsigned int impl> struct triplet long int j; T val; + triplet(long int i, long int j, T val) + { + row() = i; + col() = j; + value() = val; + } + + triplet() + { + } + long int & row() { return i; @@ -53,11 +69,38 @@ template<typename T, unsigned int impl> struct triplet * \tparam Mi implementation * */ -template<typename T,typename id_t ,typename Mi = Eigen::SparseMatrix<T,0,id_t>> +template<typename T,typename id_t ,typename Mi DEFAULT_MATRIX> class SparseMatrix { +public: + + //! Triplet implementation id + typedef boost::mpl::int_<-1> triplet_impl; + + //! Triplet type + typedef triplet<T,-1> triplet_type; + + openfpm::vector<triplet_type> stub_vt; + int stub_i; + T stub_t; + +public: + + + SparseMatrix(size_t N1, size_t N2) {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl;} + SparseMatrix() {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl;} + openfpm::vector<triplet_type> & getMatrixTriplets() {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl; return stub_vt;} + const int & getMat() const {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl; return stub_i;} + int & getMat() {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl; return stub_i;} + void resize(size_t row, size_t col) {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl;} + T operator()(id_t i, id_t j) {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl; return stub_t;} + bool save(const std::string & file) const {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl;return true;} + bool load(const std::string & file) {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl; return false;} + T getValue(size_t r, size_t c) {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl; return stub_i;} }; +#ifdef HAVE_EIGEN #include "SparseMatrix_Eigen.hpp" +#endif #endif /* OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_HPP_ */ diff --git a/src/Solvers/umfpack_solver.hpp b/src/Solvers/umfpack_solver.hpp index c7690a89eb01043340459ab2b80adf00e32efcaa..9964728cdbb900b0685184b4b4a5254f31f55b8b 100644 --- a/src/Solvers/umfpack_solver.hpp +++ b/src/Solvers/umfpack_solver.hpp @@ -8,11 +8,20 @@ #ifndef OPENFPM_NUMERICS_SRC_SOLVERS_UMFPACK_SOLVER_HPP_ #define OPENFPM_NUMERICS_SRC_SOLVERS_UMFPACK_SOLVER_HPP_ +#define UMFPACK_NONE 0 + +#define SOLVER_NOOPTION 0 +#define SOLVER_PRINT_RESIDUAL_NORM_INFINITY 1 +#define SOLVER_PRINT_DETERMINANT 2 + +#ifdef HAVE_EIGEN + +/////// Compiled with EIGEN support + #include "Vector/Vector.hpp" #include "Eigen/UmfPackSupport" #include <Eigen/SparseLU> -#define UMFPACK_NONE 0 template<typename T> class umfpack_solver @@ -25,9 +34,6 @@ public: } }; -#define SOLVER_NOOPTION 0 -#define SOLVER_PRINT_RESIDUAL_NORM_INFINITY 1 -#define SOLVER_PRINT_DETERMINANT 2 template<> class umfpack_solver<double> @@ -100,6 +106,41 @@ public: } }; +#else + +/////// Compiled without EIGEN support + +#include "Vector/Vector.hpp" + +template<typename T> +class umfpack_solver +{ +public: + + template<typename impl> static Vector<T> solve(const SparseMatrix<T,impl> & A, const Vector<T> & b) + { + std::cerr << __FILE__ << ":" << __LINE__ << " Error Umfpack only support double precision" << "/n"; + } +}; + + +template<> +class umfpack_solver<double> +{ + +public: + + template<typename impl> static Vector<double> solve(SparseMatrix<double,int,impl> & A, const Vector<double> & b, size_t opt = UMFPACK_NONE) + { + std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use umfpack you must compile OpenFPM with linear algebra support" << "/n"; + + Vector<double> x; + + return x; + } +}; + +#endif #endif /* OPENFPM_NUMERICS_SRC_SOLVERS_UMFPACK_SOLVER_HPP_ */ diff --git a/src/Vector/Vector.hpp b/src/Vector/Vector.hpp index 21e5765d21354bee983811b57caddbc8654efdd1..61973c8dc715f28c806d43055003e8d83ca3817f 100644 --- a/src/Vector/Vector.hpp +++ b/src/Vector/Vector.hpp @@ -8,8 +8,6 @@ #ifndef OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_HPP_ #define OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_HPP_ -#include <Eigen/Sparse> - /*! \brief It store one row value of a vector * * Given a row, store a value @@ -22,18 +20,50 @@ class rval }; -/* ! \brief Sparse Matrix implementation - * - * \tparam T Type of the sparse Matrix - * \tparam impl implementation +#ifdef HAVE_EIGEN +#include <Eigen/Sparse> +#define DEFAULT_VECTOR = Eigen::Matrix<T, Eigen::Dynamic, 1> +#else +#define DEFAULT_VECTOR = void +#endif + +/* ! \brief Sparse Matrix implementation stub object when OpenFPM is compiled with no linear algebra support * */ -template<typename T,typename impl = Eigen::Matrix<T, Eigen::Dynamic, 1> > +template<typename T,typename impl DEFAULT_VECTOR > class Vector { + T stub; + int stub_i; + +public: + + Vector(const Vector<T> & v) {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl;} + Vector(const Vector<T> && v) {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl;} + Vector(size_t n) {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl;} + Vector() {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl;} + void resize(size_t row) {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl;} + void insert(size_t i, T val) {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl;} + inline T & insert(size_t i) {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl;return stub;} + inline const T & insert(size_t i) const {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl; return stub;} + const T & operator()(size_t i) const {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl;return stub;} + T & operator()(size_t i) {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl;return stub;} + void scatter() {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl;} + void fromFile(std::string file) {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl;} + Vector<T> & operator=(const Vector<T> & v) {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl; return *this;} + Vector<T> & operator=(const Vector<T> && v) {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl; return *this;} + int & getVec() {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl; return stub_i;} }; +#ifdef HAVE_EIGEN #include "Vector_eigen.hpp" +#endif + +#ifdef HAVE_PETSC +#define OFPM_PETSC_VEC Vec #include "Vector_petsc.hpp" +#else +#define OFPM_PETSC_VEC void +#endif #endif /* OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_HPP_ */ diff --git a/src/Vector/Vector_eigen.hpp b/src/Vector/Vector_eigen.hpp index 00d746729bdbe141176537ac53b10741fa84990e..66c9333eea373d27e656dea2d8c3d1212c1bd6b3 100644 --- a/src/Vector/Vector_eigen.hpp +++ b/src/Vector/Vector_eigen.hpp @@ -11,12 +11,12 @@ #include <type_traits> #include "util/mul_array_extents.hpp" #include <fstream> -#include "Vector_eigen_util.hpp" #include "Grid/staggered_dist_grid.hpp" #include "Space/Ghost.hpp" #include "FiniteDifference/util/common.hpp" #include <boost/mpl/vector_c.hpp> #include <unordered_map> +#include "Vector_util.hpp" #define EIGEN_RVAL 1 diff --git a/src/Vector/Vector_eigen_util.hpp b/src/Vector/Vector_eigen_util.hpp deleted file mode 100644 index aebff63faa823cf7cbbafa30068f28b1c2b08218..0000000000000000000000000000000000000000 --- a/src/Vector/Vector_eigen_util.hpp +++ /dev/null @@ -1,124 +0,0 @@ -/* - * Vector_util.hpp - * - * Created on: Dec 7, 2015 - * Author: i-bird - */ - -#ifndef OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_EIGEN_UTIL_HPP_ -#define OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_EIGEN_UTIL_HPP_ - -#include "Grid/grid_dist_key.hpp" -#include "Space/Shape/HyperCube.hpp" - -/*! \brief Copy scalar elements - * - * \tparam copy_type Type that should be copied - * \tparam T property id to copy - * \tparam Ev Type of source the Vector - * \tparam sa dimensionality of the array 0 is a scalar - * - */ -template<typename copy_type, typename T, typename Ev, typename Eqs_sys, int sa> -struct copy_ele_sca_array -{ - template<typename Grid> inline static void copy(Grid & grid_dst, const grid_dist_key_dx<Eqs_sys::dims> & key, const Ev & x,size_t lin_id, size_t base_id, size_t gs_size) - { - grid_dst.template get<T::value>(key) = x(lin_id * Eqs_sys::nvar + base_id); - } -}; - -/*! \brief Copy 1D array elements - * - * spacialization in case of 1D array - * - * \tparam copy_type Type that should be copied - * \tparam T property id to copy - * \tparam Ev Type of source the Vector - * - */ -template<typename copy_type, typename T, typename Ev, typename Eqs_sys> -struct copy_ele_sca_array<copy_type,T,Ev,Eqs_sys,1> -{ - template<typename Grid> inline static void copy(Grid & grid_dst, const grid_dist_key_dx<Eqs_sys::dims> & key, const Ev & x,size_t lin_id, size_t base_id, size_t gs_size) - { - for (size_t i = 0 ; i < std::extent<copy_type>::value ; i++) - { - grid_dst.template get<T::value>(key)[i] = x(lin_id * Eqs_sys::nvar + base_id + i); - } - } -}; - - -/*! \brief this class is a functor for "for_each" algorithm - * - * This class is a functor for "for_each" algorithm. For each - * element of the boost::vector the operator() is called. - * Is mainly used to copy from the Vector to the grid target - * \note properties can be scalars or arrays of C++ primitives - * - * \tparam Eqs_sys System of equation information - * \tparam S type of destination grid - * \tparam Ev type of vector - * - */ -template<typename Eqs_sys, typename S, typename Ev> -struct copy_ele -{ - //! destination grid element - const grid_dist_key_dx<Eqs_sys::dims> key; - - //! destination grid - S & grid_dst; - - //! source element inside the Eigen vector - size_t lin_id; - - //! counter - size_t prp_id; - - //! It is basically the number of grid points the Eigen vector has inside - size_t gs_size; - - //! source Eigen vector - const Ev & x; - - /*! \brief constructor - * - * It define the copy parameters. - * - * \param key destination position - * \param grid_dst grid destination - * \param v Source Eigen vector - * - */ - inline copy_ele(const grid_dist_key_dx<Eqs_sys::dims> & key, S & grid_dst, const Ev & x, size_t lin_id, size_t gs_size) - :key(key),grid_dst(grid_dst),lin_id(lin_id),prp_id(0),gs_size(gs_size),x(x){}; - - -#ifdef SE_CLASS1 - /*! \brief Constructor - * - * Calling this constructor produce an error. This class store the reference of the object, - * this mean that the object passed must not be a temporal object - * - */ - inline copy_ele(grid_dist_key_dx<Eqs_sys::dims> & key, S & grid_dst, Ev && x) - :key(key),grid_dst(grid_dst),x(x) - {std::cerr << "Error: " <<__FILE__ << ":" << __LINE__ << " Passing a temporal object";}; -#endif - - //! It call the copy function for each property - template<typename T> - inline void operator()(T& t) - { - // This is the type of the object we have to copy - typedef typename boost::mpl::at_c<typename S::value_type::type,T::value>::type copy_type; - - copy_ele_sca_array<copy_type,T,Ev,Eqs_sys,std::is_array<copy_type>::value>::copy(grid_dst,key,x,lin_id,prp_id,gs_size); - prp_id += array_extents<copy_type>::mul(); - } -}; - - -#endif /* OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_EIGEN_UTIL_HPP_ */ diff --git a/src/Vector/Vector_unit_tests.hpp b/src/Vector/Vector_unit_tests.hpp index b3ee5d4bd6cba41f88c903872a0dac21f7989655..70d8cab740a189536e424a043889676365af3d66 100644 --- a/src/Vector/Vector_unit_tests.hpp +++ b/src/Vector/Vector_unit_tests.hpp @@ -139,7 +139,7 @@ BOOST_AUTO_TEST_CASE(vector_petsc_parallel) // 3 Processors 9x9 Matrix to invert - Vector<double,Vec> v(9); + Vector<double,OFPM_PETSC_VEC> v(9); if (vcl.getProcessUnitID() == 0) { @@ -161,7 +161,7 @@ BOOST_AUTO_TEST_CASE(vector_petsc_parallel) v.insert(8,8); } - Vector<double,Vec> v3; + Vector<double,OFPM_PETSC_VEC> v3; v3 = v; if (vcl.getProcessUnitID() == 0) @@ -217,7 +217,7 @@ BOOST_AUTO_TEST_CASE(vector_petsc_parallel) BOOST_REQUIRE_EQUAL(v3(2),9); } - Vec & v2 = v.getVec(); + auto & v2 = v.getVec(); } BOOST_AUTO_TEST_SUITE_END() diff --git a/src/util/grid_dist_testing.hpp b/src/util/grid_dist_testing.hpp index f49359d6114876ad9206466721ddf4216f59da33..13db51486e37af8f6517c1ed0f74c792633162c4 100644 --- a/src/util/grid_dist_testing.hpp +++ b/src/util/grid_dist_testing.hpp @@ -9,6 +9,7 @@ #define OPENFPM_NUMERICS_SRC_UTIL_GRID_DIST_TESTING_HPP_ #include "data_type/scalar.hpp" +#include "Grid/grid_dist_key.hpp" template<unsigned int dim> class grid_dist_testing diff --git a/src/util/util_num.hpp b/src/util/util_num.hpp index a51e2d42b0c962ec2caa3c501ec5760d2a3f6478..5d953e472cc80993919446eecad24a6416bb3676 100644 --- a/src/util/util_num.hpp +++ b/src/util/util_num.hpp @@ -10,6 +10,7 @@ #include "util/common.hpp" #include "grid_dist_testing.hpp" +#include "Grid/grid_dist_id.hpp" template<typename T, typename Sfinae = void> struct is_const_field: std::false_type {};