From 8868bc4095c096242ba419cd7da28ed8b9b003b4 Mon Sep 17 00:00:00 2001
From: Pietro Incardona <i-bird@private-incardon-3.mpi-cbg.de>
Date: Wed, 3 Aug 2016 17:14:28 +0200
Subject: [PATCH] Fixing numeric with latest ExtPreAlloc
---
configure.ac | 6 +
src/Equations/stoke_flow_eq_3d.hpp | 1 -
src/FiniteDifference/FDScheme.hpp | 44 +--
src/FiniteDifference/eq.hpp | 6 +
src/Makefile.am | 11 +-
src/Matrix/SparseMatrix.hpp | 5 +-
src/Matrix/SparseMatrix_Eigen.hpp | 7 +-
src/Matrix/SparseMatrix_petsc.hpp | 5 +-
src/Matrix/SparseMatrix_unit_tests.hpp | 2 +-
.../Vector/vector_dist_operators.hpp | 110 ++++++-
.../vector_dist_operators_apply_kernel.hpp | 274 ++++++++++--------
...dist_operators_apply_kernel_unit_tests.hpp | 124 +++++++-
.../vector_dist_operators_functions.hpp | 72 +++++
.../vector_dist_operators_unit_tests.hpp | 210 +++++++++++++-
src/Solvers/petsc_solver.hpp | 44 ++-
src/Solvers/umfpack_solver.hpp | 9 +-
src/Vector/Vector.hpp | 7 +-
src/Vector/Vector_petsc.hpp | 11 +-
18 files changed, 741 insertions(+), 207 deletions(-)
diff --git a/configure.ac b/configure.ac
index f070ad14..e319637a 100755
--- a/configure.ac
+++ b/configure.ac
@@ -143,6 +143,12 @@ AC_CHECK_LIB(rt, clock_gettime, [AC_DEFINE([HAVE_CLOCK_GETTIME],[],[Have clock t
OPT_LIBS="$OPT_LIBS -lrt"
])
+##########
+
+## Check for LIBHILBERT
+
+AX_LIB_HILBERT([],[echo "Cannot detect libhilbert, use the --with-libhilbert option if it is not installed in the default location"
+ exit 210])
####### include OpenFPM_numerics include path"
diff --git a/src/Equations/stoke_flow_eq_3d.hpp b/src/Equations/stoke_flow_eq_3d.hpp
index a355dc6c..886b6b38 100644
--- a/src/Equations/stoke_flow_eq_3d.hpp
+++ b/src/Equations/stoke_flow_eq_3d.hpp
@@ -12,7 +12,6 @@
constexpr unsigned int v[] = {0,1,2};
constexpr unsigned int P = 3;
-constexpr unsigned int ic = 3;
typedef Field<v[x],lid_nn_3d> v_x;
typedef Field<v[y],lid_nn_3d> v_y;
diff --git a/src/FiniteDifference/FDScheme.hpp b/src/FiniteDifference/FDScheme.hpp
index 30b0877f..7686e227 100644
--- a/src/FiniteDifference/FDScheme.hpp
+++ b/src/FiniteDifference/FDScheme.hpp
@@ -454,13 +454,32 @@ public:
*/
template<typename T> void impose(const T & op , typename Sys_eqs::stype num ,long int id ,const long int (& start)[Sys_eqs::dims], const long int (& stop)[Sys_eqs::dims], bool skip_first = false)
{
- // add padding to start and stop
- grid_key_dx<Sys_eqs::dims> start_k = grid_key_dx<Sys_eqs::dims>(start);
- grid_key_dx<Sys_eqs::dims> stop_k = grid_key_dx<Sys_eqs::dims>(stop);
+ grid_key_dx<Sys_eqs::dims> start_k;
+ grid_key_dx<Sys_eqs::dims> stop_k;
- auto it = g_map.getSubDomainIterator(start_k,stop_k);
+ bool increment = false;
+ if (skip_first == true)
+ {
+ start_k = grid_key_dx<Sys_eqs::dims>(start);
+ stop_k = grid_key_dx<Sys_eqs::dims>(start);
+
+ auto it = g_map.getSubDomainIterator(start_k,stop_k);
+
+ if (it.isNext() == true)
+ increment++;
+ }
+
+ // add padding to start and stop
+ start_k = grid_key_dx<Sys_eqs::dims>(start);
+ stop_k = grid_key_dx<Sys_eqs::dims>(stop);
+
+ auto it = g_map.getSubDomainIterator(start_k,stop_k);
+
+ if (increment == true)
+ ++it;
+
+ impose(op,num,id,it);
- impose(op,num,id,it,skip_first);
}
/*! \brief Impose an operator
@@ -478,10 +497,8 @@ public:
* \param it_d iterator that define where you want to impose
*
*/
- template<typename T> void impose(const T & op , typename Sys_eqs::stype num ,long int id ,grid_dist_iterator_sub<Sys_eqs::dims,typename g_map_type::d_grid> it_d, bool skip_first = false)
+ template<typename T> void impose(const T & op , typename Sys_eqs::stype num ,long int id ,grid_dist_iterator_sub<Sys_eqs::dims,typename g_map_type::d_grid> it_d)
{
- Vcluster & v_cl = create_vcluster();
-
openfpm::vector<triplet> & trpl = A.getMatrixTriplets();
auto it = it_d;
@@ -489,20 +506,9 @@ public:
std::unordered_map<long int,float> cols;
- bool is_first = skip_first;
-
// iterate all the grid points
while (it.isNext())
{
- if (is_first == true && v_cl.getProcessUnitID() == 0)
- {
- ++it;
- is_first = false;
- continue;
- }
- else
- is_first = false;
-
// get the position
auto key = it.get();
diff --git a/src/FiniteDifference/eq.hpp b/src/FiniteDifference/eq.hpp
index 40925175..183437ea 100644
--- a/src/FiniteDifference/eq.hpp
+++ b/src/FiniteDifference/eq.hpp
@@ -114,4 +114,10 @@ inline size_t mat_factor(size_t nvar, size_t sz, const size_t ord)
return nvar;
}
+#include "mul.hpp"
+#include "Average.hpp"
+#include "Derivative.hpp"
+#include "sum.hpp"
+#include "Laplacian.hpp"
+
#endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_EQ_HPP_ */
diff --git a/src/Makefile.am b/src/Makefile.am
index dc9b7b8f..e4f8130d 100755
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,12 +1,17 @@
-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)
+LINKLIBS = $(LIBHILBERT_LIB) $(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 = $(LIBHILBERT_INCLUDE) $(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) -lparmetis -lmetis
-nobase_include_HEADERS = PSE/Kernels.hpp PSE/Kernels_test_util.hpp Operators/Vector/vector_dist_operators_extensions.hpp Operators/Vector/vector_dist_operators.hpp Operators/Vector/vector_dist_operators_apply_kernel.hpp Operators/Vector/vector_dist_operators_functions.hpp
+nobase_include_HEADERS = Matrix/SparseMatrix.hpp Matrix/SparseMatrix_Eigen.hpp Matrix/SparseMatrix_petsc.hpp \
+Vector/Vector_eigen.hpp Vector/Vector_petsc.hpp Vector/Vector_util.hpp Vector/Vector.hpp \
+Solvers/umfpack_solver.hpp Solvers/petsc_solver.hpp \
+util/petsc_util.hpp util/linalgebra_lib.hpp util/util_num.hpp util/grid_dist_testing.hpp \
+FiniteDifference/Average.hpp FiniteDifference/Derivative.hpp FiniteDifference/eq.hpp FiniteDifference/FDScheme.hpp FiniteDifference/Laplacian.hpp FiniteDifference/mul.hpp FiniteDifference/sum.hpp FiniteDifference/util/common.hpp \
+PSE/Kernels.hpp PSE/Kernels_test_util.hpp Operators/Vector/vector_dist_operators_extensions.hpp Operators/Vector/vector_dist_operators.hpp Operators/Vector/vector_dist_operators_apply_kernel.hpp Operators/Vector/vector_dist_operators_functions.hpp Operators/Vector/vector_dist_operator_assign.hpp
.cu.o :
$(NVCC) $(NVCCFLAGS) -o $@ -c $<
diff --git a/src/Matrix/SparseMatrix.hpp b/src/Matrix/SparseMatrix.hpp
index 09d79567..0371e504 100644
--- a/src/Matrix/SparseMatrix.hpp
+++ b/src/Matrix/SparseMatrix.hpp
@@ -8,6 +8,9 @@
#ifndef OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_HPP_
#define OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_HPP_
+#include "config/config.h"
+#include "util/linalgebra_lib.hpp"
+
#ifdef HAVE_EIGEN
#include <Eigen/Sparse>
#define DEFAULT_MATRIX = EIGEN_BASE
@@ -93,7 +96,7 @@ public:
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;}
+ void resize(size_t row, size_t col, size_t row_n, size_t col_n) {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;}
diff --git a/src/Matrix/SparseMatrix_Eigen.hpp b/src/Matrix/SparseMatrix_Eigen.hpp
index c919690d..1786fd17 100644
--- a/src/Matrix/SparseMatrix_Eigen.hpp
+++ b/src/Matrix/SparseMatrix_Eigen.hpp
@@ -212,16 +212,13 @@ public:
*/
bool save(const std::string & file) const
{
- std::vector<size_t> pap_prp;
+ size_t pap_prp;
Packer<decltype(trpl),HeapMemory>::packRequest(trpl,pap_prp);
- // Calculate how much preallocated memory we need to pack all the objects
- size_t req = ExtPreAlloc<HeapMemory>::calculateMem(pap_prp);
-
// allocate the memory
HeapMemory pmem;
- pmem.allocate(req);
+ pmem.allocate(pap_prp);
ExtPreAlloc<HeapMemory> mem(pap_prp,pmem);
//Packing
diff --git a/src/Matrix/SparseMatrix_petsc.hpp b/src/Matrix/SparseMatrix_petsc.hpp
index 1a9b6b21..d129624d 100644
--- a/src/Matrix/SparseMatrix_petsc.hpp
+++ b/src/Matrix/SparseMatrix_petsc.hpp
@@ -8,7 +8,7 @@
#ifndef OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_PETSC_HPP_
#define OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_PETSC_HPP_
-
+#include "util/petsc_util.hpp"
#include "Vector/map_vector.hpp"
#include <boost/mpl/int.hpp>
#include <petscmat.h>
@@ -213,7 +213,8 @@ public:
~SparseMatrix()
{
// Destroy the matrix
- PETSC_SAFE_CALL(MatDestroy(&mat));
+ if (is_openfpm_init() == true)
+ {PETSC_SAFE_CALL(MatDestroy(&mat));}
}
/*! \brief Get the Matrix triplets bugger
diff --git a/src/Matrix/SparseMatrix_unit_tests.hpp b/src/Matrix/SparseMatrix_unit_tests.hpp
index 5a65fa41..7a9cedea 100644
--- a/src/Matrix/SparseMatrix_unit_tests.hpp
+++ b/src/Matrix/SparseMatrix_unit_tests.hpp
@@ -407,7 +407,7 @@ BOOST_AUTO_TEST_CASE(sparse_matrix_eigen_petsc_solve)
}
// Get the petsc Matrix
- Mat & matp = sm.getMat();
+ sm.getMat();
petsc_solver<double> solver;
diff --git a/src/Operators/Vector/vector_dist_operators.hpp b/src/Operators/Vector/vector_dist_operators.hpp
index a9b554b1..7adb7900 100644
--- a/src/Operators/Vector/vector_dist_operators.hpp
+++ b/src/Operators/Vector/vector_dist_operators.hpp
@@ -11,6 +11,7 @@
#include "Vector/vector_dist.hpp"
#define PROP_POS (unsigned int)-1
+#define PROP_CUSTOM (unsigned int)-2
#define VECT_SUM 1
#define VECT_SUB 2
@@ -23,12 +24,9 @@
#define VECT_APPLYKER_IN_GEN 10
#define VECT_APPLYKER_OUT_GEN 11
#define VECT_APPLYKER_REDUCE_GEN 12
-#define VECT_APPLYKER_MULTI_IN 13
-#define VECT_APPLYKER_MULTI_OUT 14
-#define VECT_APPLYKER_MULTI_REDUCE 15
-#define VECT_APPLYKER_MULTI_IN_GEN 16
-#define VECT_APPLYKER_MULTI_OUT_GEN 17
-#define VECT_APPLYKER_MULTI_REDUCE_GEN 18
+#define VECT_APPLYKER_IN_SIM 13
+#define VECT_APPLYKER_OUT_SIM 14
+#define VECT_APPLYKER_REDUCE_SIM 15
#define VECT_NORM 56
#define VECT_NORM2 57
@@ -67,6 +65,8 @@
#define VECT_PMUL 91
#define VECT_SUB_UNI 92
+#define VECT_SUM_REDUCE 93
+
/*! \brief has_init check if a type has defined a
* method called init
@@ -338,10 +338,25 @@ class vector_dist_expression
public:
+ typedef vector vtype;
+
+ //! Property id of the point
+ static const unsigned int prop = prp;
+
vector_dist_expression(vector & v)
:v(v)
{}
+ /*! \brief Return the vector on which is acting
+ *
+ * It return the vector used in getVExpr, to get this object
+ *
+ */
+ vector & getVector()
+ {
+ return v;
+ }
+
/*! \brief This function must be called before value
*
* it initialize the expression if needed
@@ -375,7 +390,7 @@ public:
{
auto key = it.get();
- v.template getProp<prp>(key) = v_exp.value(key);
+ pos_or_prop<vector,prp>::value(v,key) = v_exp.value(key);
++it;
}
@@ -383,7 +398,6 @@ public:
return v;
}
-
/*! \brief Fill the vector property with the evaluated expression
*
* \param v_exp expression to evaluate
@@ -399,7 +413,7 @@ public:
{
auto key = it.get();
- v.template getProp<prp>(key) = v_exp.value(key);
+ pos_or_prop<vector,prp>::value(v,key) = v_exp.value(key);
++it;
}
@@ -420,7 +434,7 @@ public:
{
auto key = it.get();
- v.template getProp<prp>(key) = d;
+ pos_or_prop<vector,prp>::value(v,key) = d;
++it;
}
@@ -454,7 +468,7 @@ class vector_dist_expression<prp,double>
public:
- inline vector_dist_expression(double & d)
+ inline vector_dist_expression(const double & d)
:d(d)
{}
@@ -477,6 +491,79 @@ public:
}
};
+/*! \brief Main class that encapsulate a float constant
+ *
+ * \param prp no meaning
+ *
+ */
+template<unsigned int prp>
+class vector_dist_expression<prp,float>
+{
+ float d;
+
+public:
+
+ typedef float vtype;
+
+ inline vector_dist_expression(const float & d)
+ :d(d)
+ {}
+
+ /*! \brief This function must be called before value
+ *
+ * it initialize the expression if needed
+ *
+ */
+ inline void init() const
+ {}
+
+ /*! \brief Evaluate the expression
+ *
+ * It just return the velue set in the constructor
+ *
+ */
+ inline float value(const vect_dist_key_dx & k) const
+ {
+ return d;
+ }
+};
+
+
+
+/*! \brief Main class that encapsulate a float constant
+ *
+ * \param prp no meaning
+ *
+ */
+template<typename T>
+class vector_dist_expression<PROP_CUSTOM,T>
+{
+
+public:
+
+ typedef float vtype;
+
+ inline vector_dist_expression()
+ {}
+
+ /*! \brief This function must be called before value
+ *
+ * it initialize the expression if needed
+ *
+ */
+ inline void init() const
+ {}
+
+ /*! \brief Evaluate the expression
+ *
+ * It just return the velue set in the constructor
+ *
+ */
+ inline T value(const vect_dist_key_dx & k) const
+ {
+ return T(0.0);
+ }
+};
/* \brief sum two distributed vector expression
*
@@ -1010,5 +1097,6 @@ operator/(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist
#include "vector_dist_operators_apply_kernel.hpp"
#include "vector_dist_operators_functions.hpp"
#include "vector_dist_operators_extensions.hpp"
+#include "Operators/Vector/vector_dist_operator_assign.hpp"
#endif /* OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_HPP_ */
diff --git a/src/Operators/Vector/vector_dist_operators_apply_kernel.hpp b/src/Operators/Vector/vector_dist_operators_apply_kernel.hpp
index bd10cd81..59bbc478 100644
--- a/src/Operators/Vector/vector_dist_operators_apply_kernel.hpp
+++ b/src/Operators/Vector/vector_dist_operators_apply_kernel.hpp
@@ -35,6 +35,31 @@ struct apply_kernel_rtype<exp,false>
};
+/*! \brief Meta-function to return a compatible zero-element
+ *
+ *
+ */
+template<typename rtype>
+struct set_zero
+{
+ static rtype create()
+ {
+ return 0.0;
+ }
+};
+
+template<unsigned int dim, typename T>
+struct set_zero<Point<dim,T>>
+{
+ static Point<dim,T> create()
+ {
+ Point<dim,T> ret;
+
+ ret.zero();
+ return ret;
+ }
+};
+
/*! \brief Apply the kernel to particle differently that is a number or is an expression
*
*
@@ -42,10 +67,10 @@ struct apply_kernel_rtype<exp,false>
template<typename T, typename vector, typename exp,typename NN_type, typename Kernel, typename rtype, bool is_exp=is_expression<T>::value>
struct apply_kernel_is_number_or_expression
{
- inline static rtype apply(const vector & vd, NN_type & cl, const exp & v_exp, const vect_dist_key_dx & key, Kernel & lker)
+ inline static typename std::remove_reference<rtype>::type apply(const vector & vd, NN_type & cl, const exp & v_exp, const vect_dist_key_dx & key, Kernel & lker)
{
// accumulator
- rtype pse = static_cast<rtype>(0.0);
+ typename std::remove_reference<rtype>::type pse = set_zero<typename std::remove_reference<rtype>::type>::create();
// position of particle p
Point<vector::dims,typename vector::stype> p = vd.getPos(key);
@@ -80,28 +105,21 @@ struct apply_kernel_is_number_or_expression
}
};
-template<typename T, typename vector, typename exp,typename NN_type, typename Kernel, typename rtype>
-struct apply_kernel_is_number_or_expression<T,vector,exp,NN_type,Kernel,rtype,false>
-{
- /*! \brief Apply the kernel
- *
- * \param cl Neighborhood of particles
- * \param v_exp vector expression
- * \param key particle id
- * \param lker kernel function
- *
- */
- inline static rtype apply(const vector & vd, NN_type & cl, const exp & v_exp, const vect_dist_key_dx & key, Kernel & lker)
+/*! \brief Apply the kernel to particle differently that is a number or is an expression
+ *
+ *
+ */
+template<typename vector, typename exp,typename NN_type, typename Kernel, typename rtype>
+struct apply_kernel_is_number_or_expression_sim
+{
+ inline static typename std::remove_reference<rtype>::type apply(const vector & vd, NN_type & cl, const vect_dist_key_dx & key, Kernel & lker)
{
// accumulator
- rtype pse = 0;
-
- // Get f(x) at the position of the particle p
- const rtype & prp_p = v_exp.value(key);
+ typename std::remove_reference<rtype>::type pse = set_zero<typename std::remove_reference<rtype>::type>::create();
// position of particle p
- auto & p = vd.getPos(key);
+ Point<vector::dims,typename vector::stype> p = vd.getPos(key);
// Get the neighborhood of the particle
auto NN = cl.template getNNIterator<NO_CHECK>(cl.getCell(p));
@@ -113,16 +131,10 @@ struct apply_kernel_is_number_or_expression<T,vector,exp,NN_type,Kernel,rtype,fa
// given by the Near particle, exclude itself
if (nnp != key.getKey())
{
- // Get f(x) at the position of the particle p
- const rtype & prp_q = v_exp.value(nnp);
-
- // Get f(x) at the position of the particle
- auto & q = vd.getPos(nnp);
-
- // W(x-y)
- auto ker = lker.value(p,q,prp_p,prp_q);
+ // position of the particle q
+ Point<vector::dims,typename vector::stype> q = vd.getPos(nnp);
- pse += ker;
+ pse += lker.value(p,q);
}
// Next particle
@@ -134,6 +146,7 @@ struct apply_kernel_is_number_or_expression<T,vector,exp,NN_type,Kernel,rtype,fa
};
+
/*! \brief Apply the kernel to particle differently that is a number or is an expression
*
*
@@ -141,64 +154,14 @@ struct apply_kernel_is_number_or_expression<T,vector,exp,NN_type,Kernel,rtype,fa
template<typename T, typename vector, typename exp,typename NN_type, typename Kernel, typename rtype, bool is_exp=is_expression<T>::value>
struct apply_kernel_is_number_or_expression_gen
{
- inline static rtype apply(const vector & vd, NN_type & cl, const exp & v_exp, const vect_dist_key_dx & key, Kernel & lker)
+ inline static typename std::remove_reference<rtype>::type apply(const vector & vd, NN_type & cl, const exp & v_exp, const vect_dist_key_dx & key, Kernel & lker)
{
// accumulator
- rtype pse = static_cast<rtype>(0.0);
-
- // position of particle p
- Point<vector::dims,typename vector::stype> p = vd.getPos(key);
+ typename std::remove_reference<rtype>::type pse = set_zero<typename std::remove_reference<rtype>::type>::create();
// property of the particle x
rtype prp_p = v_exp.value(key);
- // Get the neighborhood of the particle
- auto NN = cl.template getNNIterator<NO_CHECK>(cl.getCell(p));
- while(NN.isNext())
- {
- auto nnp = NN.get();
-
- // Calculate contribution given by the kernel value at position p,
- // given by the Near particle, exclude itself
- if (nnp != key.getKey())
- {
- // property of the particle x
- rtype prp_q = v_exp.value(nnp);
-
- // position of the particle q
- Point<vector::dims,typename vector::stype> q = vd.getPos(nnp);
-
- pse += lker.value(p,q,prp_p,prp_q);
- }
-
- // Next particle
- ++NN;
- }
-
- return pse;
- }
-};
-
-template<typename T, typename vector, typename exp,typename NN_type, typename Kernel, typename rtype>
-struct apply_kernel_is_number_or_expression_gen<T,vector,exp,NN_type,Kernel,rtype,false>
-{
-
- /*! \brief Apply the kernel
- *
- * \param cl Neighborhood of particles
- * \param v_exp vector expression
- * \param key particle id
- * \param lker kernel function
- *
- */
- inline static rtype apply(const vector & vd, NN_type & cl, const exp & v_exp, const vect_dist_key_dx & key, Kernel & lker)
- {
- // accumulator
- rtype pse = 0;
-
- // Get f(x) at the position of the particle p
- const rtype & prp_p = v_exp.value(key);
-
// position of particle p
auto & p = vd.getPos(key);
@@ -212,16 +175,10 @@ struct apply_kernel_is_number_or_expression_gen<T,vector,exp,NN_type,Kernel,rtyp
// given by the Near particle, exclude itself
if (nnp != key.getKey())
{
- // Get f(x) at the position of the particle p
- const rtype & prp_q = v_exp.value(nnp);
-
- // Get f(x) at the position of the particle
- auto & q = vd.getPos(nnp);
-
- // W(x-y)
- auto ker = lker.value(vd,p,q,prp_p,prp_q);
+ // property of the particle x
+ rtype prp_q = v_exp.value(nnp);
- pse += ker;
+ pse += lker.value(key.getKey(),nnp,prp_p,prp_q,vd);
}
// Next particle
@@ -232,6 +189,7 @@ struct apply_kernel_is_number_or_expression_gen<T,vector,exp,NN_type,Kernel,rtyp
}
};
+
/*! \brief Apply kernel operation
*
* \tparam exp1 expression1
@@ -273,13 +231,12 @@ public:
* \param key where to evaluate the expression
*
*/
- inline rtype value(const vect_dist_key_dx & key) const
+ inline typename std::remove_reference<rtype>::type value(const vect_dist_key_dx & key) const
{
return apply_kernel_is_number_or_expression<decltype(o1.value(key)),vector_orig,exp1,NN,Kernel,rtype>::apply(vd,cl,o1,key,ker);
}
};
-
/*! \brief Apply kernel operation
*
* \tparam exp1 expression1
@@ -287,59 +244,94 @@ public:
*
*/
template <typename exp1,typename vector_type>
-class vector_dist_expression_op<exp1,vector_type,VECT_APPLYKER_REDUCE>
+class vector_dist_expression_op<exp1,vector_type,VECT_APPLYKER_IN_SIM>
{
typedef typename boost::mpl::at<vector_type,boost::mpl::int_<0>>::type NN;
typedef typename boost::mpl::at<vector_type,boost::mpl::int_<1>>::type Kernel;
typedef typename boost::mpl::at<vector_type,boost::mpl::int_<2>>::type vector_orig;
- const exp1 o1;
NN & cl;
Kernel & ker;
const vector_orig & vd;
- typedef typename apply_kernel_rtype<decltype(o1.value(vect_dist_key_dx(0)))>::rtype rtype;
+ typedef typename apply_kernel_rtype<decltype(std::declval<Kernel>().value(Point<vector_orig::dims,typename vector_orig::stype>(0.0), Point<vector_orig::dims,typename vector_orig::stype>(0.0) ) )>::rtype rtype;
- // calculated value
- mutable rtype val;
public:
- vector_dist_expression_op(const exp1 & o1, NN & cl, Kernel & ker, const vector_orig & vd)
- :o1(o1),cl(cl),ker(ker),vd(vd)
+ /*! \brief This function must be called before value
+ *
+ * it initialize the expression if needed
+ *
+ */
+ inline void init() const
+ {
+ }
+
+ vector_dist_expression_op(NN & cl, Kernel & ker, const vector_orig & vd)
+ :cl(cl),ker(ker),vd(vd)
{}
- /*! \brief Initialize the calculation
+ /*! \brief Evaluate the expression
*
+ * \param key where to evaluate the expression
*
*/
- void init() const
+ inline typename std::remove_reference<rtype>::type value(const vect_dist_key_dx & key) const
{
- val = 0.0;
+ return apply_kernel_is_number_or_expression_sim<vector_orig,exp1,NN,Kernel,rtype>::apply(vd,cl,key,ker);
+ }
+};
+
+
+/*! \brief Apply kernel operation
+ *
+ * \tparam exp1 expression1
+ * \tparam NN list
+ *
+ */
+template <typename exp1,typename vector_type>
+class vector_dist_expression_op<exp1,vector_type,VECT_APPLYKER_IN_GEN>
+{
+ typedef typename boost::mpl::at<vector_type,boost::mpl::int_<0>>::type NN;
+ typedef typename boost::mpl::at<vector_type,boost::mpl::int_<1>>::type Kernel;
+ typedef typename boost::mpl::at<vector_type,boost::mpl::int_<2>>::type vector_orig;
- auto it = vd.getDomainIterator();
+ const exp1 o1;
+ NN & cl;
+ Kernel & ker;
+ const vector_orig & vd;
- while (it.isNext())
- {
- auto key = it.get();
+ typedef typename apply_kernel_rtype<decltype(o1.value(vect_dist_key_dx(0)))>::rtype rtype;
- val += apply_kernel_is_number_or_expression<decltype(o1.value(key)),vector_orig,exp1,NN,Kernel,rtype>::apply(vd,cl,o1,key,ker);
+public:
- ++it;
- }
+ /*! \brief This function must be called before value
+ *
+ * it initialize the expression if needed
+ *
+ */
+ inline void init() const
+ {
+ o1.init();
}
+ vector_dist_expression_op(const exp1 & o1, NN & cl, Kernel & ker, const vector_orig & vd)
+ :o1(o1),cl(cl),ker(ker),vd(vd)
+ {}
+
/*! \brief Evaluate the expression
*
* \param key where to evaluate the expression
*
*/
- inline rtype value(const vect_dist_key_dx & key) const
+ inline typename std::remove_reference<rtype>::type value(const vect_dist_key_dx & key) const
{
- return val;
+ return apply_kernel_is_number_or_expression_gen<decltype(o1.value(key)),vector_orig,exp1,NN,Kernel,rtype>::apply(vd,cl,o1,key,ker);
}
};
+
///////////////////////////////////////////// Apply kernel operator ////
////////////////////////////////////////////////////////////////////////
@@ -360,6 +352,7 @@ applyKernel_in(const vector_dist_expression_op<exp1,exp2,op1> & va, vector_type
return exp_sum;
}
+
/* \brief Apply kernel expression
*
* \param va vector expression one
@@ -369,14 +362,69 @@ applyKernel_in(const vector_dist_expression_op<exp1,exp2,op1> & va, vector_type
*
*/
template<typename exp1, typename exp2, unsigned int op1, typename NN, typename Kernel, typename vector_type>
-inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,boost::mpl::vector<NN,Kernel,vector_type>,VECT_APPLYKER_REDUCE>
-applyKernel_reduce(const vector_dist_expression_op<exp1,exp2,op1> & va, vector_type & vd, NN & cl, Kernel & ker)
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,boost::mpl::vector<NN,Kernel,vector_type>,VECT_APPLYKER_IN_GEN>
+applyKernel_in_gen(const vector_dist_expression_op<exp1,exp2,op1> & va, vector_type & vd, NN & cl, Kernel & ker)
{
- vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,boost::mpl::vector<NN,Kernel,vector_type>,VECT_APPLYKER_REDUCE> exp_sum(va,cl,ker,vd);
+ vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,boost::mpl::vector<NN,Kernel,vector_type>,VECT_APPLYKER_IN_GEN> exp_sum(va,cl,ker,vd);
return exp_sum;
}
+//////////////////////////////////////// For vector expression ///////////////////////
+
+/* \brief Apply kernel expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<unsigned int prp, typename NN, typename Kernel, typename vector_type>
+inline vector_dist_expression_op<vector_dist_expression<prp,vector_type>,boost::mpl::vector<NN,Kernel,vector_type>,VECT_APPLYKER_IN>
+applyKernel_in(const vector_dist_expression<prp,vector_type> & va, vector_type & vd, NN & cl, Kernel & ker)
+{
+ vector_dist_expression_op<vector_dist_expression<prp,vector_type>,boost::mpl::vector<NN,Kernel,vector_type>,VECT_APPLYKER_IN> exp_sum(va,cl,ker,vd);
+
+ return exp_sum;
+}
+
+
+/* \brief Apply kernel expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<unsigned int prp, typename NN, typename Kernel, typename vector_type>
+inline vector_dist_expression_op<vector_dist_expression<prp,vector_type>,boost::mpl::vector<NN,Kernel,vector_type>,VECT_APPLYKER_IN_GEN>
+applyKernel_in_gen(const vector_dist_expression<prp,vector_type> & va, vector_type & vd, NN & cl, Kernel & ker)
+{
+ vector_dist_expression_op<vector_dist_expression<prp,vector_type>,boost::mpl::vector<NN,Kernel,vector_type>,VECT_APPLYKER_IN_GEN> exp_sum(va,cl,ker,vd);
+
+ return exp_sum;
+}
+
+
+/* \brief Apply kernel expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<typename NN, typename Kernel, typename vector_type>
+inline vector_dist_expression_op<void,boost::mpl::vector<NN,Kernel,vector_type>,VECT_APPLYKER_IN_SIM>
+applyKernel_in_sim(vector_type & vd, NN & cl, Kernel & ker)
+{
+ vector_dist_expression_op<void,boost::mpl::vector<NN,Kernel,vector_type>,VECT_APPLYKER_IN_SIM> exp_sum(cl,ker,vd);
+
+ return exp_sum;
+}
+
#endif /* OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_APPLY_KERNEL_HPP_ */
diff --git a/src/Operators/Vector/vector_dist_operators_apply_kernel_unit_tests.hpp b/src/Operators/Vector/vector_dist_operators_apply_kernel_unit_tests.hpp
index ff345fab..efd7ab9c 100644
--- a/src/Operators/Vector/vector_dist_operators_apply_kernel_unit_tests.hpp
+++ b/src/Operators/Vector/vector_dist_operators_apply_kernel_unit_tests.hpp
@@ -51,6 +51,9 @@ template <typename vector,typename Kernel, typename NN_type> bool check_values_a
base2 += vd.template getProp<C>(p);
+ if (base1 != base2)
+ std::cout << "CAZZO " << base1 << " " << base2 << std::endl;
+
ret &= base1 == base2;
++it2;
@@ -176,6 +179,58 @@ template <typename vector,typename Kernel, typename NN_type> bool check_values_a
return ret;
}
+template <typename vector,typename Kernel, typename NN_type> bool check_values_apply_kernel3(vector & vd, Kernel & ker, NN_type & NN)
+{
+ bool ret = true;
+ auto it = vd.getDomainIterator();
+
+ // ### WIKI 13 ###
+ //
+ // Check that apply kernel work
+ //
+ auto it2 = vd.getDomainIterator();
+ while (it2.isNext())
+ {
+ Point<3,float> base2 = 0.0;
+ auto p = it2.get();
+
+ Point<3,float> xp = vd.getPos(p);
+
+ Point<3,float> base1 = vd.template getProp<VA>(p);
+
+ Point<3,float> prp_x = vd.template getProp<VC>(p);
+
+ // For each neighborhood particle
+ auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.getPos(p)));
+
+ while (Np.isNext())
+ {
+ // Neighborhood particle q
+ auto q = Np.get();
+
+ if (q == p.getKey()) {++Np; continue;};
+
+ // position q
+ Point<3,float> xq = vd.getPos(q);
+
+ Point<3,float> prp_y = vd.template getProp<VC>(q);
+
+ base2 += ker.value(xp,xq,prp_x,prp_y);
+
+ ++Np;
+ }
+
+ base2 += vd.template getProp<VC>(p);
+
+ ret &= base1 == base2;
+
+ ++it2;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+
+ return ret;
+}
template <typename vector,typename Kernel, typename NN_type> bool check_values_apply_kernel2_reduce(vector & vd, Kernel & ker, NN_type & NN)
{
@@ -238,6 +293,51 @@ template <typename vector,typename Kernel, typename NN_type> bool check_values_a
return ret;
}
+template <typename vector,typename Kernel, typename NN_type> bool check_values_apply_kernel3_reduce(vector & vd, Kernel & ker, NN_type & NN, const Point<2,float> & p)
+{
+ bool ret = true;
+ auto it = vd.getDomainIterator();
+
+ Point<2,float> base2 = 0.0;
+
+ auto it2 = vd.getDomainIterator();
+ while (it2.isNext())
+ {
+ auto p = it2.get();
+
+ Point<3,float> xp = vd.getPos(p);
+
+ Point<2,float> ker_accu = 0.0;
+
+ // For each neighborhood particle
+ auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.getPos(p)));
+
+ while (Np.isNext())
+ {
+ // Neighborhood particle q
+ auto q = Np.get();
+
+ if (q == p.getKey()) {++Np; continue;};
+
+ // position q
+ Point<3,float> xq = vd.getPos(q);
+
+ ker_accu += ker.value(xp,xq);
+
+ ++Np;
+ }
+
+ base2 += ker_accu;
+
+ ++it2;
+ }
+
+ BOOST_REQUIRE_EQUAL(p.get(0),base2.get(0));
+ BOOST_REQUIRE_EQUAL(p.get(1),base2.get(1));
+
+ return ret;
+}
+
BOOST_AUTO_TEST_CASE( vector_dist_operators_apply_kernel_test )
{
if (create_vcluster().getProcessingUnits() > 3)
@@ -257,7 +357,6 @@ BOOST_AUTO_TEST_CASE( vector_dist_operators_apply_kernel_test )
vector_dist<3,float,aggregate<float,float,float,VectorS<3,float>,VectorS<3,float>,VectorS<3,float>,float>> vd(512,box,bc,ghost);
auto vA = getV<A>(vd);
- auto vB = getV<B>(vd);
auto vC = getV<C>(vd);
auto vVA = getV<VA>(vd);
@@ -281,23 +380,34 @@ BOOST_AUTO_TEST_CASE( vector_dist_operators_apply_kernel_test )
vVA = applyKernel_in(2.0*vVC + vVB ,vd,cl,ker) + vVC;
check_values_apply_kernel2(vd,ker,cl);
- vA = applyKernel_reduce(vVC * vVB + norm(vVB),vd,cl,ker) + vC;
+ vA = rsum(applyKernel_in(vVC * vVB + norm(vVB),vd,cl,ker),vd) + vC;
check_values_apply_kernel_reduce(vd,ker,cl);
- vVA = applyKernel_reduce(2.0*vVC + vVB ,vd,cl,ker) + vVC;
+ vVA = rsum(applyKernel_in(2.0*vVC + vVB ,vd,cl,ker),vd) + vVC;
check_values_apply_kernel2_reduce(vd,ker,cl);
-/* vA = applyKernel_in_gen(vVC * vVB + norm(vVB),vd,cl,ker) + vC;
+ vA = applyKernel_in_gen(vVC * vVB + norm(vVB),vd,cl,ker) + vC;
check_values_apply_kernel(vd,ker,cl);
vVA = applyKernel_in_gen(2.0*vVC + vVB ,vd,cl,ker) + vVC;
check_values_apply_kernel2(vd,ker,cl);
- vA = applyKernel_reduce_gen(vVC * vVB + norm(vVB),vd,cl,ker) + vC;
+ vA = rsum(applyKernel_in_gen(vVC * vVB + norm(vVB),vd,cl,ker),vd) + vC;
check_values_apply_kernel_reduce(vd,ker,cl);
- vVA = applyKernel_reduce_gen(2.0*vVC + vVB ,vd,cl,ker) + vVC;
- check_values_apply_kernel2_reduce(vd,ker,cl);*/
+ vVA = rsum(applyKernel_in_gen(2.0*vVC + vVB ,vd,cl,ker),vd) + vVC;
+ check_values_apply_kernel2_reduce(vd,ker,cl);
+
+ // Check it compile the code is the same
+ vVA = applyKernel_in_gen(vVC,vd,cl,ker) + vVC;
+ check_values_apply_kernel3(vd,ker,cl);
+
+ vVA = applyKernel_in (vVC,vd,cl,ker) + vVC;
+ check_values_apply_kernel3(vd,ker,cl);
+
+ // just check it compile
+ Point<2,float> p = rsum(applyKernel_in_sim(vd,cl,ker),vd).get();
+ check_values_apply_kernel3_reduce(vd,ker,cl,p);
}
diff --git a/src/Operators/Vector/vector_dist_operators_functions.hpp b/src/Operators/Vector/vector_dist_operators_functions.hpp
index c53919ba..c0d9dfb3 100644
--- a/src/Operators/Vector/vector_dist_operators_functions.hpp
+++ b/src/Operators/Vector/vector_dist_operators_functions.hpp
@@ -212,4 +212,76 @@ fun_name(double d, const vector_dist_expression_op<exp1,exp2,op1> & va)\
CREATE_VDIST_ARG2_FUNC(pmul,pmul,VECT_PMUL)
+
+////////// Special function reduce /////////////////////////
+
+
+template <typename exp1, typename vector_type>
+class vector_dist_expression_op<exp1,vector_type,VECT_SUM_REDUCE>
+{
+ const exp1 o1;
+
+ typedef typename apply_kernel_rtype<decltype(o1.value(vect_dist_key_dx(0)))>::rtype rtype;
+
+ // calculated value
+ mutable typename std::remove_reference<rtype>::type val;
+
+ const vector_type & vd;
+
+public:
+
+ vector_dist_expression_op(const exp1 & o1, const vector_type & vd)
+ :o1(o1), vd(vd)
+ {}
+
+ inline void init() const
+ {
+ o1.init();
+
+ val = 0.0;
+
+ auto it = vd.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ val += o1.value(key);
+
+ ++it;
+ }
+ }
+
+ inline typename std::remove_reference<rtype>::type get()
+ {
+ init();
+ return value(vect_dist_key_dx(0));
+ }
+
+ template<typename r_type= typename std::remove_reference<rtype>::type > inline r_type value(const vect_dist_key_dx & key) const
+ {
+ return val;
+ }
+};
+
+
+template<typename exp1, typename exp2_, unsigned int op1, typename vector_type>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2_,op1>,vector_type,VECT_SUM_REDUCE>
+rsum(const vector_dist_expression_op<exp1,exp2_,op1> & va, const vector_type & vd)
+{
+ vector_dist_expression_op<vector_dist_expression_op<exp1,exp2_,op1>,vector_type,VECT_SUM_REDUCE> exp_sum(va,vd);
+
+ return exp_sum;
+}
+
+template<unsigned int prp1, typename v1, typename vector_type>
+inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_type,VECT_SUM_REDUCE>
+rsum(const vector_dist_expression<prp1,v1> & va, const vector_type & vd)
+{
+ vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_type,VECT_SUM_REDUCE> exp_sum(va,vd);
+
+ return exp_sum;
+}
+
+
#endif /* OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_FUNCTIONS_HPP_ */
diff --git a/src/Operators/Vector/vector_dist_operators_unit_tests.hpp b/src/Operators/Vector/vector_dist_operators_unit_tests.hpp
index 208f7c83..7a0e9040 100644
--- a/src/Operators/Vector/vector_dist_operators_unit_tests.hpp
+++ b/src/Operators/Vector/vector_dist_operators_unit_tests.hpp
@@ -646,7 +646,8 @@ template <typename vector> void fill_values(vector & v)
}
}
-typedef vector_dist<3,float,aggregate<float,float,float,VectorS<3,float>,VectorS<3,float>,VectorS<3,float>>> vector_type;
+
+typedef vector_dist<3,float,aggregate<float,float,float,VectorS<3,float>,VectorS<3,float>,VectorS<3,float>,float>> vector_type;
//! Exponential kernel
struct exp_kernel
@@ -671,17 +672,37 @@ struct exp_kernel
return (pA + pB) * exp(dist * dist / var);
}
- inline Point<3,float> value(vector_type & vd1, vector_type & vd2, size_t p, size_t q, float pA, float pB)
+ inline float value(size_t p, size_t q, float pA, float pB, const vector_type & vd1)
{
- float dist = norm(p-q);
+ Point<3,float> pp = vd1.getPos(p);
+ Point<3,float> pq = vd1.getPos(q);
+
+ float dist = norm(pp-pq);
return (pA + pB) * exp(dist * dist / var);
}
+
+ inline Point<3,float> value(size_t p, size_t q, const Point<3,float> & pA, const Point<3,float> & pB , const vector_type & vd1)
+ {
+ Point<3,float> pp = vd1.getPos(p);
+ Point<3,float> pq = vd1.getPos(q);
+
+ float dist = norm(pp-pq);
+
+ return (pA + pB) * exp(dist * dist / var);
+ }
+
+ inline Point<2,float> value(const Point<3,float> & p, const Point<3,float> & q)
+ {
+ float dist = norm(p-q);
+
+ return exp(dist * dist / var);
+ }
};
////////////////////////////////////////////////////////////////////////////
-BOOST_AUTO_TEST_SUITE( vector_dist_test )
+BOOST_AUTO_TEST_SUITE( vector_dist_operators_test )
BOOST_AUTO_TEST_CASE( vector_dist_operators_test )
@@ -882,12 +903,12 @@ BOOST_AUTO_TEST_CASE( vector_dist_operators_test )
vA = norm(vPOS);
vVA = vPOS + 2.0;
vVA = 2.0 + vPOS;
-// vVA = vPOS + vPOS;
+ vVA = vPOS + vPOS;
vVA = vPOS - 2.0f;
vVA = 2.0 - vPOS;
-// vVA = vPOS - vPOS;
+ vVA = vPOS - vPOS;
-/* vVA = vPOS * 2.0;
+ vVA = vPOS * 2.0;
vVA = 2.0 * vPOS;
vVA = vPOS * vPOS;
@@ -914,11 +935,184 @@ BOOST_AUTO_TEST_CASE( vector_dist_operators_test )
vVA = vPOS / (vPOS + vPOS);
vVA = (vPOS + vPOS) / vPOS;
- vVA = (vPOS + vPOS) / (vPOS + vPOS);*/
+ vVA = (vPOS + vPOS) / (vPOS + vPOS);
}
#include "vector_dist_operators_apply_kernel_unit_tests.hpp"
+void reset(vector_dist<3,float,aggregate<float,float,float,float,float,float,float,float>> & v1)
+{
+ auto it = v1.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto p = it.get();
+
+ v1.template getProp<0>(p) = -1.0;
+ v1.template getProp<1>(p) = -1.0;
+ v1.template getProp<2>(p) = -1.0;
+ v1.template getProp<3>(p) = -1.0;
+ v1.template getProp<4>(p) = -1.0;
+ v1.template getProp<5>(p) = -1.0;
+ v1.template getProp<6>(p) = -1.0;
+ v1.template getProp<7>(p) = -1.0;
+
+ ++it;
+ }
+}
+
+template<unsigned int i> void loop_check(vector_dist<3,float,aggregate<float,float,float,float,float,float,float,float>> & v1, float check)
+{
+ auto it = v1.getDomainIterator();
+ bool ret = true;
+
+ while (it.isNext())
+ {
+ auto p = it.get();
+
+ ret &= v1.template getProp<i>(p) == check;
+
+ ++it;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+}
+
+void check(vector_dist<3,float,aggregate<float,float,float,float,float,float,float,float>> & v1, size_t i)
+{
+ float check = -1.0;
+
+ if (0 < i) check = 0.0;
+ else check = -1.0;
+ loop_check<0>(v1,check);
+
+ if (1 < i) check = 1.0;
+ else check = -1.0;
+ loop_check<1>(v1,check);
+
+ if (2 < i) check = 2.0;
+ else check = -1.0;
+ loop_check<2>(v1,check);
+
+ if (3 < i) check = 3.0;
+ else check = -1.0;
+ loop_check<3>(v1,check);
+
+ if (4 < i) check = 4.0;
+ else check = -1.0;
+ loop_check<4>(v1,check);
+
+ if (5 < i) check = 5.0;
+ else check = -1.0;
+ loop_check<5>(v1,check);
+
+ if (6 < i) check = 6.0;
+ else check = -1.0;
+ loop_check<6>(v1,check);
+
+ if (7 < i) check = 7.0;
+ else check = -1.0;
+ loop_check<7>(v1,check);
+}
+
+BOOST_AUTO_TEST_CASE( vector_dist_operators_assign_test )
+{
+ if (create_vcluster().getProcessingUnits() > 3)
+ return;
+
+ Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
+
+ // Boundary conditions
+ size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
+
+ // ghost
+ Ghost<3,float> ghost(0.05);
+
+ // vector type
+ typedef vector_dist<3,float,aggregate<float,float,float,float,float,float,float,float>> vtype;
+
+ vector_dist<3,float,aggregate<float,float,float,float,float,float,float,float>> vd(100,box,bc,ghost);
+
+ auto v1 = getV<0>(vd);
+ auto v2 = getV<1>(vd);
+ auto v3 = getV<2>(vd);
+ auto v4 = getV<3>(vd);
+ auto v5 = getV<4>(vd);
+ auto v6 = getV<5>(vd);
+ auto v7 = getV<6>(vd);
+ auto v8 = getV<7>(vd);
+
+ reset(vd);
+
+ assign(v1,0.0,
+ v2,1.0);
+
+ check(vd,2);
+
+ reset(vd);
+
+ assign(v1,0.0,
+ v2,1.0,
+ v3,2.0);
+
+ check(vd,3);
+
+ reset(vd);
+
+ assign(v1,0.0,
+ v2,1.0,
+ v3,2.0,
+ v4,3.0);
+
+ check(vd,4);
+
+ reset(vd);
+
+ assign(v1,0.0,
+ v2,1.0,
+ v3,2.0,
+ v4,3.0,
+ v5,4.0);
+
+ check(vd,5);
+
+ reset(vd);
+
+ assign(v1,0.0,
+ v2,1.0,
+ v3,2.0,
+ v4,3.0,
+ v5,4.0,
+ v6,5.0);
+
+ check(vd,6);
+
+ reset(vd);
+
+ assign(v1,0.0,
+ v2,1.0,
+ v3,2.0,
+ v4,3.0,
+ v5,4.0,
+ v6,5.0,
+ v7,6.0);
+
+ check(vd,7);
+
+ reset(vd);
+
+ assign(v1,0.0,
+ v2,1.0,
+ v3,2.0,
+ v4,3.0,
+ v5,4.0,
+ v6,5.0,
+ v7,6.0,
+ v8,7.0);
+
+ check(vd,8);
+}
+
BOOST_AUTO_TEST_SUITE_END()
diff --git a/src/Solvers/petsc_solver.hpp b/src/Solvers/petsc_solver.hpp
index 04d0c91e..2a82cb2c 100644
--- a/src/Solvers/petsc_solver.hpp
+++ b/src/Solvers/petsc_solver.hpp
@@ -96,15 +96,6 @@ class petsc_solver<double>
openfpm::vector<itError> res;
};
- // KSP Relative tolerance
-// PetscReal rtol;
-
- // KSP Absolute tolerance
-// PetscReal abstol;
-
- // KSP dtol tolerance to determine that the method diverge
-// PetscReal dtol;
-
// KSP Maximum number of iterations
PetscInt maxits;
@@ -227,7 +218,7 @@ class petsc_solver<double>
yn.add(bench.get(i).smethod);
}
- size_t n_int = 300;
+ size_t n_int = maxits;
// calculate dt
@@ -303,18 +294,22 @@ class petsc_solver<double>
itError erri(err);
erri.it = it;
- //
- size_t old_size = pts->bench.last().res.size();
- pts->bench.last().res.resize(it+1);
+ itError err_fill;
- if (old_size > 0)
- {
- for (size_t i = old_size ; i < it-1 ; i++)
- pts->bench.last().res.get(i) = pts->bench.last().res.get(i-1);
- }
+ size_t old_size = pts->bench.last().res.size();
+ pts->bench.last().res.resize(it+1);
+
+ if (old_size > 0)
+ err_fill = pts->bench.last().res.get(old_size-1);
+ else
+ err_fill = erri;
+
+ for (long int i = old_size ; i < (long int)it ; i++)
+ pts->bench.last().res.get(i) = err_fill;
+
+ // Add the error per iteration
+ pts->bench.last().res.get(it) = erri;
- // Add the error per iteration
- pts->bench.last().res.add(erri);
pts->progress(it);
@@ -683,7 +678,6 @@ class petsc_solver<double>
void initKSP()
{
PETSC_SAFE_CALL(KSPCreate(PETSC_COMM_WORLD,&ksp));
-// PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
}
/*! \brief initialize the KSP object for solver testing
@@ -693,10 +687,8 @@ class petsc_solver<double>
void initKSPForTest()
{
PETSC_SAFE_CALL(KSPCreate(PETSC_COMM_WORLD,&ksp));
-// PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
-// PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol,abstol,dtol,300));
- setMaxIter(300);
+ setMaxIter(maxits);
// Disable convergence check
PETSC_SAFE_CALL(KSPSetConvergenceTest(ksp,KSPConvergedSkip,NULL,NULL));
@@ -718,6 +710,7 @@ public:
}
petsc_solver()
+ :maxits(300)
{
initKSP();
@@ -725,7 +718,7 @@ public:
solvs.add(std::string(KSPBCGS));
// solvs.add(std::string(KSPIBCGS)); <--- Produce invalid argument
- solvs.add(std::string(KSPFBCGS));
+// solvs.add(std::string(KSPFBCGS));
// solvs.add(std::string(KSPFBCGSR)); <--- Nan production problems
solvs.add(std::string(KSPBCGSL));
solvs.add(std::string(KSPGMRES));
@@ -831,6 +824,7 @@ public:
void setMaxIter(PetscInt n)
{
PetscOptionsSetValue("-ksp_max_it",std::to_string(n).c_str());
+ maxits = n;
}
/*! For the BiCGStab(L) it define the number of search directions
diff --git a/src/Solvers/umfpack_solver.hpp b/src/Solvers/umfpack_solver.hpp
index e1b1f4e1..8a05f0ce 100644
--- a/src/Solvers/umfpack_solver.hpp
+++ b/src/Solvers/umfpack_solver.hpp
@@ -28,9 +28,9 @@ class umfpack_solver
{
public:
- template<typename impl> static Vector<T> solve(const SparseMatrix<T,impl> & A, const Vector<T> & b)
+ template<unsigned int impl, typename id_type> static Vector<T> solve(const SparseMatrix<T,id_type,impl> & A, const Vector<T> & b)
{
- std::cerr << "Error Umfpack only suppor double precision" << "/n";
+ std::cerr << "Error Umfpack only support double precision, and int ad id type" << "/n";
}
};
@@ -84,6 +84,9 @@ public:
{
// decomposition failed
std::cout << __FILE__ << ":" << __LINE__ << " solver failed" << "\n";
+
+ x.scatter();
+
return x;
}
@@ -136,7 +139,7 @@ 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)
+ template<unsigned int impl, typename id_type> static Vector<double> solve(SparseMatrix<double,id_type,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";
diff --git a/src/Vector/Vector.hpp b/src/Vector/Vector.hpp
index b11b131e..24ef44f6 100644
--- a/src/Vector/Vector.hpp
+++ b/src/Vector/Vector.hpp
@@ -8,6 +8,7 @@
#ifndef OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_HPP_
#define OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_HPP_
+#include "config/config.h"
#include "util/linalgebra_lib.hpp"
/*! \brief It store one row value of a vector
@@ -43,11 +44,11 @@ class Vector
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(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;}
Vector(size_t n, size_t n_row_local) {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 resize(size_t row, size_t row_n) {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;}
@@ -56,7 +57,7 @@ public:
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;}
+ Vector<T> & operator=(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;}
};
diff --git a/src/Vector/Vector_petsc.hpp b/src/Vector/Vector_petsc.hpp
index d3b20416..b24981b8 100644
--- a/src/Vector/Vector_petsc.hpp
+++ b/src/Vector/Vector_petsc.hpp
@@ -131,18 +131,19 @@ public:
* \param v vector to copy
*
*/
- Vector(const Vector<T,PETSC_BASE> && v)
+ Vector(Vector<T,PETSC_BASE> && v)
{
this->operator=(v);
}
- /*! \brief Destoroy the vector
+ /*! \brief Destroy the vector
*
*
*/
~Vector()
{
- PETSC_SAFE_CALL(VecDestroy(&v));
+ if (is_openfpm_init() == true)
+ {PETSC_SAFE_CALL(VecDestroy(&v));}
}
/*! \brief Create a vector with n elements
@@ -330,7 +331,7 @@ public:
// Fill the index and construct the map
size_t k = 0;
- for (size_t i = low ; i < high ; i++)
+ for (size_t i = low ; i < (size_t)high ; i++)
{
row_val.template get<row_id>(k) = i;
map[i] = k;
@@ -358,7 +359,7 @@ public:
* \param v vector to copy
*
*/
- Vector<T,PETSC_BASE> & operator=(const Vector<T,PETSC_BASE> && v)
+ Vector<T,PETSC_BASE> & operator=(Vector<T,PETSC_BASE> && v)
{
map.swap(v.map);
row_val.swap(v.row_val);
--
GitLab