diff --git a/src/FiniteDifference/FDScheme.hpp b/src/FiniteDifference/FDScheme.hpp
index 0488cec1cccae752b8e4b24951f0a00f17a78bda..64d9ad79fd50fabf7e79210292f92fc202ccbfc1 100644
--- a/src/FiniteDifference/FDScheme.hpp
+++ b/src/FiniteDifference/FDScheme.hpp
@@ -365,7 +365,9 @@ public:
// sync the ghost
+ g_map.write("g_map_before_ghost.vtk");
g_map.template ghost_get<0>();
+ g_map.write("g_map_after_ghost.vtk");
// Create a CellDecomposer and calculate the spacing
@@ -512,7 +514,7 @@ public:
// copy the vector
for (size_t i = 0; i < row_b; i++)
- B.get(i) = b.get(i);
+ B.insert(i,b.get(i));
return B;
}
diff --git a/src/FiniteDifference/eq.hpp b/src/FiniteDifference/eq.hpp
index 6d85cc9e6dbff7d37666ac9ffad86b05df36221c..f4af2d96074dca9ffa1628647a55fd97e7ffeb47 100644
--- a/src/FiniteDifference/eq.hpp
+++ b/src/FiniteDifference/eq.hpp
@@ -11,9 +11,6 @@
#define EQS_FIELD 0
#define EQS_POS 1
-#define STAGGERED_GRID 1
-#define NORMAL_GRID 0
-
//#define PERIODIC true
//#define NON_PERIODIC false
diff --git a/src/FiniteDifference/eq_unit_test.hpp b/src/FiniteDifference/eq_unit_test.hpp
index 842692ea3dc09ef5dc66fbeef4d14adad51ead21..6ebf2ab87d16538e09dc3c9d1bbed524c05dacb5 100644
--- a/src/FiniteDifference/eq_unit_test.hpp
+++ b/src/FiniteDifference/eq_unit_test.hpp
@@ -142,13 +142,13 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
//! [lid-driven cavity 2D]
// Domain, a rectangle
- Box<2,float> domain({0.0,0.0},{3.0,1.0});
+ Box<2,float> domain({0.0,0.0},{1.0,1.0});
// Ghost (Not important in this case but required)
Ghost<2,float> g(0.01);
// Grid points on x=256 and y=64
- long int sz[] = {256,64};
+ long int sz[] = {8,8};
size_t szu[2];
szu[0] = (size_t)sz[0];
szu[1] = (size_t)sz[1];
diff --git a/src/FiniteDifference/util/common.hpp b/src/FiniteDifference/util/common.hpp
index 4b383d6b795f36c720d7cb17c0a5d579eedd1b34..dc71a5e2e76509966db2b2ec0cb831a851b41965 100644
--- a/src/FiniteDifference/util/common.hpp
+++ b/src/FiniteDifference/util/common.hpp
@@ -8,6 +8,8 @@
#ifndef OPENFPM_NUMERICS_SRC_UTIL_COMMON_HPP_
#define OPENFPM_NUMERICS_SRC_UTIL_COMMON_HPP_
+#define STAGGERED_GRID 1
+#define NORMAL_GRID 0
template<typename T, typename Sfinae = void>
struct has_grid_type: std::false_type {};
diff --git a/src/Matrix/SparseMatrix_Eigen.hpp b/src/Matrix/SparseMatrix_Eigen.hpp
index 706c818cb496e0e54106903438833003819d5c40..85bb274cb3256129531908a55b43b3877a77c3a8 100644
--- a/src/Matrix/SparseMatrix_Eigen.hpp
+++ b/src/Matrix/SparseMatrix_Eigen.hpp
@@ -21,30 +21,49 @@
*
*/
template<typename T>
-struct triplet<T,EIGEN_TRIPLET>
+class triplet<T,EIGEN_TRIPLET>
{
- Eigen::Triplet<T,long int> triplet;
+ Eigen::Triplet<T,long int> triplet_ei;
+
+public:
long int & row()
{
- size_t ptr = (size_t)&triplet.row();
+ size_t ptr = (size_t)&triplet_ei.row();
return (long int &)*(long int *)ptr;
}
long int & col()
{
- size_t ptr = (size_t)&triplet.col();
+ size_t ptr = (size_t)&triplet_ei.col();
return (long int &)*(long int *)ptr;
}
T & value()
{
- size_t ptr = (size_t)&triplet.value();
+ size_t ptr = (size_t)&triplet_ei.value();
return (T &)*(T *)ptr;
}
+
+ /*! \brief Constructor from row, colum and value
+ *
+ * \param i row
+ * \param j colum
+ * \param val value
+ *
+ */
+ triplet(long int i, long int j, T val)
+ {
+ row() = i;
+ col() = j;
+ value() = val;
+ }
+
+ // Default constructor
+ triplet() {};
};
/* ! \brief Sparse Matrix implementation, that map over Eigen
@@ -73,8 +92,6 @@ private:
/*! \brief Assemble the matrix
*
- * \param trpl Matrix in form of triplets
- * \param mat Matrix to assemble
*
*/
void assemble()
@@ -187,6 +204,16 @@ public:
{
mat.resize(row,col);
}
+
+ /*! \brief Get the row i and the colum j of the Matrix
+ *
+ * \return the value of the matrix at row i colum j
+ *
+ */
+ T operator()(id_t i, id_t j)
+ {
+ return mat.coeffRef(i,j);
+ }
};
diff --git a/src/Solvers/umfpack_solver.hpp b/src/Solvers/umfpack_solver.hpp
index 826343cd036b6c1727349af46eb647c73ad5ec4e..2d5f67d30d4291851416bd6aeb8a6b09cd477029 100644
--- a/src/Solvers/umfpack_solver.hpp
+++ b/src/Solvers/umfpack_solver.hpp
@@ -46,16 +46,27 @@ public:
*/
template<typename impl> static Vector<double> solve(SparseMatrix<double,int,impl> & A, const Vector<double> & b, size_t opt = UMFPACK_NONE)
{
- Vcluster & v_cl = *global_v_cluster;
+ Vcluster & vcl = *global_v_cluster;
Vector<double> x;
+
// only master processor solve
+ Eigen::UmfPackLU<Eigen::SparseMatrix<double,0,int> > solver;
- if (v_cl.getProcessUnitID() == 0)
- {
- Eigen::UmfPackLU<Eigen::SparseMatrix<double,0,int> > solver;
+ // Collect the matrix on master
+ auto mat_ei = A.getMat();
+
+ Eigen::Matrix<double, Eigen::Dynamic, 1> x_ei;
+
+ // Collect the vector on master
+ auto b_ei = b.getVec();
- solver.compute(A.getMat());
+ // Copy b into x, this also copy the information on how to scatter back the information on x
+ x = b;
+
+ if (vcl.getProcessUnitID() == 0)
+ {
+ solver.compute(mat_ei);
if(solver.info()!=Eigen::Success)
{
@@ -64,14 +75,14 @@ public:
return x;
}
- x.getVec() = solver.solve(b.getVec());
+ x_ei = solver.solve(b_ei);
if (opt & SOLVER_PRINT_RESIDUAL_NORM_INFINITY)
{
- Vector<double> res;
- res.getVec() = A.getMat() * x.getVec() - b.getVec();
+ Eigen::Matrix<double, Eigen::Dynamic, 1> res;
+ res = mat_ei * x_ei - b_ei;
- std::cout << "Infinity norm: " << res.getVec().lpNorm<Eigen::Infinity>() << "\n";
+ std::cout << "Infinity norm: " << res.lpNorm<Eigen::Infinity>() << "\n";
}
if (opt & SOLVER_PRINT_DETERMINANT)
@@ -79,9 +90,12 @@ public:
std::cout << " Determinant: " << solver.determinant() << "\n";
}
- // Vector is only on master, scatter back the information
- x.scatter();
+ x = x_ei;
}
+
+ // Vector is only on master, scatter back the information
+ x.scatter();
+
return x;
}
};
diff --git a/src/Vector/Vector.hpp b/src/Vector/Vector.hpp
index 8277376953f5284a128bed93d79f68a59d12e055..43469dba28c9163c2aab486164ecd50beab0beaa 100644
--- a/src/Vector/Vector.hpp
+++ b/src/Vector/Vector.hpp
@@ -8,6 +8,19 @@
#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
+ *
+ *
+ */
+template<typename T, unsigned int impl>
+class rval
+{
+};
+
/* ! \brief Sparse Matrix implementation
*
diff --git a/src/Vector/Vector_eigen.hpp b/src/Vector/Vector_eigen.hpp
index f0d5c2b8cc6d8a1eeab61a895af9ea00de094d17..4ef615b6c7eb9c825515ef6ea3523038d0f11078 100644
--- a/src/Vector/Vector_eigen.hpp
+++ b/src/Vector/Vector_eigen.hpp
@@ -13,13 +13,80 @@
#include <fstream>
#include "Vector_eigen_util.hpp"
#include "Grid/staggered_dist_grid_util.hpp"
+#include "Space/Ghost.hpp"
+#include "FiniteDifference/util/common.hpp"
#include <boost/mpl/vector_c.hpp>
+#include <unordered_map>
+#define EIGEN_RVAL 1
+
+/*! \brief It store one row value of a vector
+ *
+ * Given a row, store a value
+ *
+ *
+ */
+template<typename T>
+class rval<T,EIGEN_RVAL>
+{
+ // row
+ long int r;
+
+ // value
+ T val;
+
+public:
+
+ // Get the row
+ long int & row()
+ {
+ return r;
+ }
+
+ // Get the value
+ T & value()
+ {
+ return val;
+ }
+
+ /*! \brief Default constructor
+ *
+ */
+ rval() {}
+
+ /*! \brief Constructor from row, colum and value
+ *
+ * \param i row
+ * \param val value
+ *
+ */
+ rval(long int i, T val)
+ {
+ row() = i;
+ value() = val;
+ }
+};
template<typename T>
class Vector<T,Eigen::Matrix<T, Eigen::Dynamic, 1>>
{
- Eigen::Matrix<T, Eigen::Dynamic, 1> v;
+ mutable Eigen::Matrix<T, Eigen::Dynamic, 1> v;
+
+ // row value vector
+ mutable openfpm::vector<rval<T,EIGEN_RVAL>> row_val;
+ mutable openfpm::vector<rval<T,EIGEN_RVAL>> row_val_recv;
+
+ // global to local map
+ mutable std::unordered_map<size_t,size_t> map;
+
+ // invalid
+ T invalid;
+
+ // Processors from where we gather
+ mutable openfpm::vector<size_t> prc;
+
+ //size of each chunk
+ mutable openfpm::vector<size_t> sz;
/*! \brief Check that the size of the iterators match
*
@@ -169,8 +236,90 @@ class Vector<T,Eigen::Matrix<T, Eigen::Dynamic, 1>>
}
}
+
+ /*! \brief Here we collect the full matrix on master
+ *
+ */
+ void collect() const
+ {
+ Vcluster & vcl = *global_v_cluster;
+
+ row_val_recv.clear();
+
+ // here we collect all the triplet in one array on the root node
+ vcl.SGather(row_val,row_val_recv,prc,sz,0);
+
+ if (vcl.getProcessUnitID() != 0)
+ row_val.resize(0);
+ else
+ row_val.swap(row_val_recv);
+
+ build_map();
+ }
+
+ /*! \brief Set the Eigen internal vector
+ *
+ *
+ */
+ void setEigen() const
+ {
+ // set the vector
+
+ for (size_t i = 0 ; i < row_val.size() ; i++)
+ v[row_val.get(i).row()] = row_val.get(i).value();
+ }
+
+ /*! \brief Build the map
+ *
+ *
+ */
+ void build_map() const
+ {
+ map.clear();
+
+ for (size_t i = 0 ; i < row_val.size() ; i++)
+ map[row_val.get(i).row()] = i;
+ }
+
public:
+ /*! \brief Copy the vector
+ *
+ * \param v vector to copy
+ *
+ */
+ Vector(const Vector<T> & v)
+ {
+ this->operator=(v);
+ }
+
+ /*! \brief Copy the vector
+ *
+ * \param v vector to copy
+ *
+ */
+ Vector(const Vector<T> && v)
+ {
+ this->operator=(v);
+ }
+
+ /*! \brief Create a vector with n elements
+ *
+ * \param n number of elements in the vector
+ *
+ */
+ Vector(size_t n)
+ {
+ resize(n);
+ }
+
+ /*! \brief Create a vector with 0 elements
+ *
+ */
+ Vector()
+ {
+ }
+
/*! \brief Resize the Vector
*
* \param row numbers of row
@@ -188,9 +337,44 @@ public:
* \return reference to the element vector
*
*/
- T & get(size_t i)
+ void insert(size_t i, T val)
{
- return v(i);
+ row_val.add();
+
+ // Map
+ map[i] = row_val.size()-1;
+
+ row_val.last().row() = i;
+ row_val.last().value() = val;
+ }
+
+ /*! \brief Return a reference to the vector element
+ *
+ * \warning The element must exist
+ *
+ * \param i element
+ *
+ * \return reference to the element vector
+ *
+ */
+ T & operator()(size_t i)
+ {
+ // Search if exist
+
+ std::unordered_map<size_t,size_t>::iterator it = map.find(i);
+
+ if ( it == map.end() )
+ {
+ // Does not exist
+
+ std::cerr << __FILE__ << ":" << __LINE__ << " value does not exist " << std::endl;
+
+ return invalid;
+ }
+ else
+ return row_val.get(it->second).value();
+
+ return invalid;
}
/*! \brief Get the Eigen Vector object
@@ -200,6 +384,9 @@ public:
*/
const Eigen::Matrix<T, Eigen::Dynamic, 1> & getVec() const
{
+ collect();
+ setEigen();
+
return v;
}
@@ -210,6 +397,9 @@ public:
*/
Eigen::Matrix<T, Eigen::Dynamic, 1> & getVec()
{
+ collect();
+ setEigen();
+
return v;
}
@@ -228,30 +418,28 @@ public:
}
}
- /*! \brief Collect the vector into one processor
- *
- * Eigen does not have a real parallel vector, so in order to work we have to collect
- * the vector in one processor
- *
- * This function collect the information from all the other processor into one
- *
- */
- void collect()
- {
-
- }
-
/*! \brief Scatter the vector information to the other processors
*
* Eigen does not have a real parallel vector, so in order to work we have to scatter
* the vector from one processor to the other
*
- * This function scatter the information to all the other processors
*
*/
void scatter()
{
+ row_val_recv.clear();
+ Vcluster & vcl = *global_v_cluster;
+
+ vcl.SScatter(row_val,row_val_recv,prc,sz,0);
+ // if we do not receive anything a previous collect has not been performed
+ // and so nothing is scattered
+ if (row_val_recv.size() != 0)
+ {
+ row_val.clear();
+ row_val.add(row_val_recv);
+ build_map();
+ }
}
/*! \brief Load from file
@@ -270,6 +458,51 @@ public:
inputf.close();
}
+
+ /*! \brief Copy the vector
+ *
+ * \param v vector to copy
+ *
+ */
+ Vector<T> & operator=(const Vector<T> & v)
+ {
+ prc = v.prc;
+ sz = v.sz;
+ map = v.map;
+ row_val = v.row_val;
+
+ return *this;
+ }
+
+ /*! \brief Copy the vector
+ *
+ * \param v vector to copy
+ *
+ */
+ Vector<T> & operator=(const Vector<T> && v)
+ {
+ prc = v.prc;
+ sz = v.sz;
+ map = v.map;
+ row_val = v.row_val;
+
+ return *this;
+ }
+
+ /*! \brief Copy the vector (it is used for special purpose)
+ *
+ * \warning v MUST contain at least all the elements of the vector
+ *
+ * \param v base eigen vector to copy
+ *
+ */
+ Vector<T> & operator=(Eigen::Matrix<T, Eigen::Dynamic, 1> & v)
+ {
+ for (size_t i = 0 ; i < row_val.size() ; i++)
+ row_val.get(i).value() = v(row_val.get(i).row());
+
+ return *this;
+ }
};
diff --git a/src/Vector/Vector_eigen_util.hpp b/src/Vector/Vector_eigen_util.hpp
index da16b2a7868c1b429e6fdb7c4cf08b5664c79a42..221a8600be5e5e4bb819c62ff85b9c2bb294b67e 100644
--- a/src/Vector/Vector_eigen_util.hpp
+++ b/src/Vector/Vector_eigen_util.hpp
@@ -8,6 +8,8 @@
#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
*
diff --git a/src/main.cpp b/src/main.cpp
index 57fb816dcd55bd6455b17767d16396a66b4c6868..8a63fcf886529ad092a5e275a6a60affdd764194 100755
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -4,6 +4,8 @@
#include "unit_test_init_cleanup.hpp"
+#include "Vector/Vector_unit_tests.hpp"
+#include "Matrix/SparseMatrix_unit_tests.hpp"
#include "Equations/eq_unit_tests.hpp"
#include "FiniteDifference/FDScheme_unit_tests.hpp"
#include "FiniteDifference/util/common_test.hpp"