Skip to content
Snippets Groups Projects
Commit 2fc10c99 authored by Pietro Incardona's avatar Pietro Incardona
Browse files

Project compile even without linear algebra

parent 6d7ddcff
No related branches found
No related tags found
No related merge requests found
...@@ -18,7 +18,7 @@ ...@@ -18,7 +18,7 @@
#include "Grid/comb.hpp" #include "Grid/comb.hpp"
#include "FiniteDifference/util/common.hpp" #include "FiniteDifference/util/common.hpp"
#include "util/util_num.hpp" #include "util/util_num.hpp"
#include <unordered_map>
/*! \brief Derivative second order on h (spacing) /*! \brief Derivative second order on h (spacing)
* *
......
...@@ -15,6 +15,9 @@ ...@@ -15,6 +15,9 @@
#include "data_type/scalar.hpp" #include "data_type/scalar.hpp"
#include "NN/CellList/CellDecomposer.hpp" #include "NN/CellList/CellDecomposer.hpp"
#include "Grid/staggered_dist_grid_util.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 /*! \brief Finite Differences
* *
......
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 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_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_CFLAGS = $(CUDA_CFLAGS)
numerics_LDADD = $(LINKLIBS) numerics_LDADD = $(LINKLIBS) -lmetis -lparmetis
nobase_include_HEADERS = PSE/Kernels.hpp PSE/Kernels_test_util.hpp nobase_include_HEADERS = PSE/Kernels.hpp PSE/Kernels_test_util.hpp
.cu.o : .cu.o :
......
...@@ -8,7 +8,12 @@ ...@@ -8,7 +8,12 @@
#ifndef OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_HPP_ #ifndef OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_HPP_
#define OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_HPP_ #define OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_HPP_
#ifdef HAVE_EIGEN
#include <Eigen/Sparse> #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 /*! \brief It store the non zero elements of the matrix
* *
...@@ -30,6 +35,17 @@ template<typename T, unsigned int impl> struct triplet ...@@ -30,6 +35,17 @@ template<typename T, unsigned int impl> struct triplet
long int j; long int j;
T val; T val;
triplet(long int i, long int j, T val)
{
row() = i;
col() = j;
value() = val;
}
triplet()
{
}
long int & row() long int & row()
{ {
return i; return i;
...@@ -53,11 +69,38 @@ template<typename T, unsigned int impl> struct triplet ...@@ -53,11 +69,38 @@ template<typename T, unsigned int impl> struct triplet
* \tparam Mi implementation * \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 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" #include "SparseMatrix_Eigen.hpp"
#endif
#endif /* OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_HPP_ */ #endif /* OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_HPP_ */
...@@ -8,11 +8,20 @@ ...@@ -8,11 +8,20 @@
#ifndef OPENFPM_NUMERICS_SRC_SOLVERS_UMFPACK_SOLVER_HPP_ #ifndef OPENFPM_NUMERICS_SRC_SOLVERS_UMFPACK_SOLVER_HPP_
#define 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 "Vector/Vector.hpp"
#include "Eigen/UmfPackSupport" #include "Eigen/UmfPackSupport"
#include <Eigen/SparseLU> #include <Eigen/SparseLU>
#define UMFPACK_NONE 0
template<typename T> template<typename T>
class umfpack_solver class umfpack_solver
...@@ -25,9 +34,6 @@ public: ...@@ -25,9 +34,6 @@ public:
} }
}; };
#define SOLVER_NOOPTION 0
#define SOLVER_PRINT_RESIDUAL_NORM_INFINITY 1
#define SOLVER_PRINT_DETERMINANT 2
template<> template<>
class umfpack_solver<double> class umfpack_solver<double>
...@@ -100,6 +106,41 @@ public: ...@@ -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_ */ #endif /* OPENFPM_NUMERICS_SRC_SOLVERS_UMFPACK_SOLVER_HPP_ */
...@@ -8,8 +8,6 @@ ...@@ -8,8 +8,6 @@
#ifndef OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_HPP_ #ifndef OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_HPP_
#define 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 /*! \brief It store one row value of a vector
* *
* Given a row, store a value * Given a row, store a value
...@@ -22,18 +20,50 @@ class rval ...@@ -22,18 +20,50 @@ class rval
}; };
/* ! \brief Sparse Matrix implementation #ifdef HAVE_EIGEN
* #include <Eigen/Sparse>
* \tparam T Type of the sparse Matrix #define DEFAULT_VECTOR = Eigen::Matrix<T, Eigen::Dynamic, 1>
* \tparam impl implementation #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 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" #include "Vector_eigen.hpp"
#endif
#ifdef HAVE_PETSC
#define OFPM_PETSC_VEC Vec
#include "Vector_petsc.hpp" #include "Vector_petsc.hpp"
#else
#define OFPM_PETSC_VEC void
#endif
#endif /* OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_HPP_ */ #endif /* OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_HPP_ */
...@@ -11,12 +11,12 @@ ...@@ -11,12 +11,12 @@
#include <type_traits> #include <type_traits>
#include "util/mul_array_extents.hpp" #include "util/mul_array_extents.hpp"
#include <fstream> #include <fstream>
#include "Vector_eigen_util.hpp"
#include "Grid/staggered_dist_grid.hpp" #include "Grid/staggered_dist_grid.hpp"
#include "Space/Ghost.hpp" #include "Space/Ghost.hpp"
#include "FiniteDifference/util/common.hpp" #include "FiniteDifference/util/common.hpp"
#include <boost/mpl/vector_c.hpp> #include <boost/mpl/vector_c.hpp>
#include <unordered_map> #include <unordered_map>
#include "Vector_util.hpp"
#define EIGEN_RVAL 1 #define EIGEN_RVAL 1
......
/*
* 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_ */
...@@ -139,7 +139,7 @@ BOOST_AUTO_TEST_CASE(vector_petsc_parallel) ...@@ -139,7 +139,7 @@ BOOST_AUTO_TEST_CASE(vector_petsc_parallel)
// 3 Processors 9x9 Matrix to invert // 3 Processors 9x9 Matrix to invert
Vector<double,Vec> v(9); Vector<double,OFPM_PETSC_VEC> v(9);
if (vcl.getProcessUnitID() == 0) if (vcl.getProcessUnitID() == 0)
{ {
...@@ -161,7 +161,7 @@ BOOST_AUTO_TEST_CASE(vector_petsc_parallel) ...@@ -161,7 +161,7 @@ BOOST_AUTO_TEST_CASE(vector_petsc_parallel)
v.insert(8,8); v.insert(8,8);
} }
Vector<double,Vec> v3; Vector<double,OFPM_PETSC_VEC> v3;
v3 = v; v3 = v;
if (vcl.getProcessUnitID() == 0) if (vcl.getProcessUnitID() == 0)
...@@ -217,7 +217,7 @@ BOOST_AUTO_TEST_CASE(vector_petsc_parallel) ...@@ -217,7 +217,7 @@ BOOST_AUTO_TEST_CASE(vector_petsc_parallel)
BOOST_REQUIRE_EQUAL(v3(2),9); BOOST_REQUIRE_EQUAL(v3(2),9);
} }
Vec & v2 = v.getVec(); auto & v2 = v.getVec();
} }
BOOST_AUTO_TEST_SUITE_END() BOOST_AUTO_TEST_SUITE_END()
......
...@@ -9,6 +9,7 @@ ...@@ -9,6 +9,7 @@
#define OPENFPM_NUMERICS_SRC_UTIL_GRID_DIST_TESTING_HPP_ #define OPENFPM_NUMERICS_SRC_UTIL_GRID_DIST_TESTING_HPP_
#include "data_type/scalar.hpp" #include "data_type/scalar.hpp"
#include "Grid/grid_dist_key.hpp"
template<unsigned int dim> template<unsigned int dim>
class grid_dist_testing class grid_dist_testing
......
...@@ -10,6 +10,7 @@ ...@@ -10,6 +10,7 @@
#include "util/common.hpp" #include "util/common.hpp"
#include "grid_dist_testing.hpp" #include "grid_dist_testing.hpp"
#include "Grid/grid_dist_id.hpp"
template<typename T, typename Sfinae = void> template<typename T, typename Sfinae = void>
struct is_const_field: std::false_type {}; struct is_const_field: std::false_type {};
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment