Commit 3d8ac690 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

corrections in vector_operations

parent e5d7c85d
...@@ -243,8 +243,8 @@ protected: ...@@ -243,8 +243,8 @@ protected:
bool value1set = false; bool value1set = false;
if (asmMatrix) { if (asmMatrix) {
typedef traits::range_generator<tag::row, Matrix>::type c_type; typedef mtl::traits::range_generator<tag::row, Matrix>::type c_type;
typedef traits::range_generator<tag::nz, c_type>::type ic_type; typedef mtl::traits::range_generator<tag::nz, c_type>::type ic_type;
for (size_t col = 0; col < getNumComponents(); col++) { for (size_t col = 0; col < getNumComponents(); col++) {
TEST_EXIT(getSystemMatrix(row_, col) != NULL || col != col_) TEST_EXIT(getSystemMatrix(row_, col) != NULL || col != col_)
...@@ -255,9 +255,9 @@ protected: ...@@ -255,9 +255,9 @@ protected:
// set Dirichlet-row in matrix // set Dirichlet-row in matrix
Matrix &m = getSystemMatrix(row_, col)->getBaseMatrix(); Matrix &m = getSystemMatrix(row_, col)->getBaseMatrix();
traits::row<Matrix>::type r(m); mtl::traits::row<Matrix>::type r(m);
traits::col<Matrix>::type c(m); mtl::traits::col<Matrix>::type c(m);
traits::value<Matrix>::type v(m); mtl::traits::value<Matrix>::type v(m);
c_type cursor(begin<tag::row>(m)+idx_); c_type cursor(begin<tag::row>(m)+idx_);
for (ic_type icursor(begin<tag::nz>(cursor)), icend(end<tag::nz>(cursor)); icursor != icend; ++icursor) { for (ic_type icursor(begin<tag::nz>(cursor)), icend(end<tag::nz>(cursor)); icursor != icend; ++icursor) {
...@@ -287,8 +287,8 @@ protected: ...@@ -287,8 +287,8 @@ protected:
using namespace mtl; using namespace mtl;
typedef DOFMatrix::base_matrix_type Matrix; typedef DOFMatrix::base_matrix_type Matrix;
typedef traits::range_generator<tag::row, Matrix>::type c_type; typedef mtl::traits::range_generator<tag::row, Matrix>::type c_type;
typedef traits::range_generator<tag::nz, c_type>::type ic_type; typedef mtl::traits::range_generator<tag::nz, c_type>::type ic_type;
for (size_t col = 0; col < getNumComponents(); col++) { for (size_t col = 0; col < getNumComponents(); col++) {
TEST_EXIT(getSystemMatrix(row, col) != NULL) TEST_EXIT(getSystemMatrix(row, col) != NULL)
...@@ -296,9 +296,9 @@ protected: ...@@ -296,9 +296,9 @@ protected:
Matrix &m = getSystemMatrix(row, col)->getBaseMatrix(); Matrix &m = getSystemMatrix(row, col)->getBaseMatrix();
traits::row<Matrix>::type r(m); mtl::traits::row<Matrix>::type r(m);
traits::col<Matrix>::type c(m); mtl::traits::col<Matrix>::type c(m);
traits::value<Matrix>::type v(m); mtl::traits::value<Matrix>::type v(m);
std::vector<std::vector<std::pair<DegreeOfFreedom, double> > > row_values; std::vector<std::vector<std::pair<DegreeOfFreedom, double> > > row_values;
row_values.resize(indices.size()); row_values.resize(indices.size());
...@@ -349,8 +349,8 @@ protected: ...@@ -349,8 +349,8 @@ protected:
using namespace mtl; using namespace mtl;
typedef DOFMatrix::base_matrix_type Matrix; typedef DOFMatrix::base_matrix_type Matrix;
typedef traits::range_generator<tag::row, Matrix>::type c_type; typedef mtl::traits::range_generator<tag::row, Matrix>::type c_type;
typedef traits::range_generator<tag::nz, c_type>::type ic_type; typedef mtl::traits::range_generator<tag::nz, c_type>::type ic_type;
TEST_EXIT(row_idx.size() == coefficients.size() && row_idx.size() == rhs.size() && rhs.size()>0) TEST_EXIT(row_idx.size() == coefficients.size() && row_idx.size() == rhs.size() && rhs.size()>0)
("rhs_idx, coefficients and rhs must have the same size and size >! 0\n"); ("rhs_idx, coefficients and rhs must have the same size and size >! 0\n");
...@@ -363,9 +363,9 @@ protected: ...@@ -363,9 +363,9 @@ protected:
Matrix &m = getSystemMatrix(row, col)->getBaseMatrix(); Matrix &m = getSystemMatrix(row, col)->getBaseMatrix();
traits::row<Matrix>::type r(m); mtl::traits::row<Matrix>::type r(m);
traits::col<Matrix>::type c(m); mtl::traits::col<Matrix>::type c(m);
traits::value<Matrix>::type v(m); mtl::traits::value<Matrix>::type v(m);
// erase the rows for all row-indices and set rhs values // erase the rows for all row-indices and set rhs values
for (size_t i = 0; i < coefficients.size(); i++) { for (size_t i = 0; i < coefficients.size(); i++) {
......
...@@ -9,7 +9,6 @@ ...@@ -9,7 +9,6 @@
using namespace std; using namespace std;
using namespace AMDiS; using namespace AMDiS;
using namespace vector_operations;
/** /**
* A collection of signed-dist function describing several gemoetric objects: * A collection of signed-dist function describing several gemoetric objects:
...@@ -33,6 +32,8 @@ static double signedDist2D(const WorldVector<double> x, const WorldVector<double ...@@ -33,6 +32,8 @@ static double signedDist2D(const WorldVector<double> x, const WorldVector<double
AbstractFunction<double, double> *radius, double eps) AbstractFunction<double, double> *radius, double eps)
{ {
FUNCNAME("signedDist2D"); FUNCNAME("signedDist2D");
using vector_operations::epsNorm;
WorldVector<double> x_trans; WorldVector<double> x_trans;
double norm_xy; double norm_xy;
double alpha; double alpha;
...@@ -59,6 +60,8 @@ static double signedDist3D(const WorldVector<double> x, const WorldVector<double ...@@ -59,6 +60,8 @@ static double signedDist3D(const WorldVector<double> x, const WorldVector<double
BinaryAbstractFunction<double, double, double> *radius, double eps) BinaryAbstractFunction<double, double, double> *radius, double eps)
{ {
FUNCNAME("signedDist3D"); FUNCNAME("signedDist3D");
using vector_operations::epsNorm;
WorldVector<double> x_trans; WorldVector<double> x_trans;
double norm_xyz, norm_xy; double norm_xyz, norm_xy;
double alpha, beta; double alpha, beta;
......
/** \file ValueTypes.h */
#ifndef VALUE_TYPES_H
#define VALUE_TYPES_H
#include "AMDiS.h"
#include "boost/type_traits.hpp"
using namespace AMDiS;
template<typename T> struct DOFView;
/// Type-Traits for value_type of Data-Structures
/// ________________________________________________________________________________________________
template<typename T, typename enable=void> struct ValueType { typedef T type; };
template<typename T> struct ValueType<DOFVector<T> > { typedef T type; };
template<typename T> struct ValueType<std::vector<T> > { typedef T type; };
template<typename T> struct ValueType<std::list<T> > { typedef T type; };
template<typename T> struct ValueType<std::set<T> > { typedef T type; };
template<typename T> struct ValueType<mtl::dense_vector<T> > { typedef T type; };
template<typename Derived> struct ValueType<Derived, typename boost::enable_if<boost::is_base_of<DOFView<double>, Derived > >::type> { typedef double type; };
template<typename Derived> struct ValueType<Derived, typename boost::enable_if<boost::is_base_of<DOFView<WorldVector<double> >, Derived > >::type> { typedef WorldVector<double> type; };
template<typename Derived> struct ValueType<Derived, typename boost::enable_if<boost::is_base_of<AbstractFunction<double, WorldVector<double> >, Derived > >::type> { typedef double type; };
template<typename Derived> struct ValueType<Derived, typename boost::enable_if<boost::is_base_of<AbstractFunction<WorldVector<double>, WorldVector<double> >, Derived > >::type> { typedef WorldVector<double> type; };
/// random-accessible vectors
/// ________________________________________________________________________________________________
template<typename Vector> struct is_vector : public boost::false_type {};
template<typename T> struct is_vector< std::vector<T> > : public boost::true_type {};
template<typename T> struct is_vector< mtl::dense_vector<T> > : public boost::true_type {};
template<typename T> struct is_vector< WorldVector<T> > : public boost::true_type {};
/// random-accessible matrices
/// ________________________________________________________________________________________________
template<typename Matrix> struct is_matrix : public boost::false_type {};
template<typename T, typename Param> struct is_matrix< mtl::matrix::base_matrix<T, Param> > : public boost::true_type {};
template<typename T> struct is_matrix< mtl::matrix::dense2D<T> > : public boost::true_type {};
template<typename T> struct is_matrix< mtl::matrix::compressed2D<T> > : public boost::true_type {};
template<typename T> struct is_matrix< WorldMatrix<T> > : public boost::true_type {};
#endif // VALUE_TYPES_H
...@@ -4,6 +4,7 @@ ...@@ -4,6 +4,7 @@
#define VECTOR_OPERATIONS_H #define VECTOR_OPERATIONS_H
#include "AMDiS.h" #include "AMDiS.h"
#include "ValueTypes.h"
#if HAVE_PARALLEL_DOMAIN_AMDIS #if HAVE_PARALLEL_DOMAIN_AMDIS
#include "parallel/StdMpi.h" #include "parallel/StdMpi.h"
...@@ -16,6 +17,7 @@ ...@@ -16,6 +17,7 @@
#include <string> #include <string>
#include <boost/numeric/mtl/mtl.hpp> #include <boost/numeric/mtl/mtl.hpp>
#include <boost/numeric/itl/itl.hpp> #include <boost/numeric/itl/itl.hpp>
#include <boost/utility/enable_if.hpp>
using namespace AMDiS; using namespace AMDiS;
...@@ -76,67 +78,176 @@ public: ...@@ -76,67 +78,176 @@ public:
namespace vector_operations { namespace vector_operations {
// num_rows for vector types namespace traits {
template<typename T>
size_t num_rows(WorldVector<T> &v) {
return static_cast<size_t>(v.getSize());
};
template<typename T>
size_t num_rows(const WorldVector<T> &v) {
return static_cast<size_t>(v.getSize());
};
template<typename T>
size_t num_rows(std::vector<T> &v) {
return static_cast<size_t>(v.size());
};
template<typename T>
size_t num_rows(const std::vector<T> &v) {
return static_cast<size_t>(v.size());
};
template<typename T>
size_t num_rows(T &v) {
return 1;
};
template<typename T>
size_t num_rows(const T &v) {
return 1;
};
// num_rows/num_cols for matrix types /// General declaration, used to disable unsupported types
template<typename T> template <typename Collection, class Enable = void>
size_t num_rows(WorldMatrix<T> &m) { struct num_rows {};
return static_cast<size_t>(m.getNumRows());
}; /// size implementation for STL vectors
template<typename T> template <typename Value>
size_t num_rows(const WorldMatrix<T> &m) { struct num_rows< std::vector<Value> >
return static_cast<size_t>(m.getNumRows()); {
}; typedef std::size_t type;
template<typename T> type operator()(const std::vector<Value>& v) { return v.size(); }
size_t num_cols(WorldMatrix<T> &m) { };
return static_cast<size_t>(m.getNumCols());
}; /// size implementation for (1D) arrays interpreted as vectors
template<typename T> template <typename Value, unsigned Size>
size_t num_cols(const WorldMatrix<T> &m) { struct num_rows<Value[Size]>
return static_cast<size_t>(m.getNumCols()); {
}; typedef std::size_t type;
type operator()(const Value[Size]) { return Size; }
// common interface to resize vectors };
template<typename T>
void resize(WorldVector<T> &v, size_t dim) { /// size implementation for (2D and higher) arrays interpreted as matrices
TEST_EXIT(dim==v.getSize())("WorldVectors can not be resized to the given dimension!\n"); template <typename Value, unsigned Rows, unsigned Cols>
}; struct num_rows<Value[Rows][Cols]>
template<typename T> {
void resize(std::vector<T> &v, size_t dim) { typedef std::size_t type;
v.resize(dim); type operator()(const Value[Rows][Cols]) { return Rows; }
}; };
template<typename T>
void resize(mtl::dense_vector<T> &v, size_t dim) { template <typename Value>
v.change_dim(dim); struct num_rows< WorldVector<Value> >
}; {
typedef std::size_t type;
type operator()(const WorldVector<Value>& v) { return static_cast<size_t>(v.getSize()); }
};
template <typename Value>
struct num_rows< WorldMatrix<Value> >
{
typedef std::size_t type;
type operator()(const WorldMatrix<Value>& v) { return static_cast<size_t>(v.getNumRows()); }
};
template <typename Value>
struct num_rows< Value, typename boost::enable_if< boost::is_arithmetic<Value> >::type >
{
typedef std::size_t type;
type operator()(const Value& v) { return 1; }
};
//________________________________________________________________________________
/// General declaration, used to disable unsupported types
template <typename Collection, class Enable = void>
struct num_cols {};
template <typename Value>
struct num_cols< WorldMatrix<Value> >
{
typedef std::size_t type;
type operator()(const WorldMatrix<Value>& v) { return static_cast<size_t>(v.getNumCols()); }
};
//________________________________________________________________________________
/// General declaration, used to disable unsupported types
template <typename Collection, class Enable = void>
struct resize {};
template <typename Value>
struct resize< WorldVector<Value> >
{
typedef void vector_void_type;
void operator()(WorldVector<Value>& v, size_t r) {
TEST_EXIT(Global::getGeo(WORLD) == r)
("WorldVectors can not be resized!\n");
}
};
template <typename Value>
struct resize< std::vector<Value> >
{
typedef void vector_void_type;
void operator()(std::vector<Value>& v, size_t r) {
v.resize(r);
}
};
template <typename Value>
struct resize< mtl::dense_vector<Value> >
{
typedef void vector_void_type;
void operator()(mtl::dense_vector<Value>& v, size_t r) {
v.change_dim(r);
}
};
// _________________________________________________________________________________
template <typename Value>
struct resize< WorldMatrix<Value> >
{
typedef void matrix_void_type;
void operator()(WorldMatrix<Value>& v, size_t r, size_t c) {
TEST_EXIT(Global::getGeo(WORLD) == r && Global::getGeo(WORLD) == c)
("WorldMatrices can not be resized!\n");
}
};
template <typename Value, typename Param>
struct resize< mtl::matrix::base_matrix<Value, Param> >
{
typedef void matrix_void_type;
void operator()(mtl::matrix::base_matrix<Value, Param>& v, size_t r, size_t c) {
v.change_dim(r,c);
}
};
template <typename Value>
struct resize< mtl::matrix::dense2D<Value> >
{
typedef void matrix_void_type;
void operator()(mtl::matrix::dense2D<Value>& v, size_t r, size_t c) {
v.change_dim(r,c);
}
};
template <typename Value>
struct resize< mtl::matrix::compressed2D<Value> >
{
typedef void matrix_void_type;
void operator()(mtl::matrix::compressed2D<Value>& v, size_t r, size_t c) {
v.change_dim(r,c);
}
};
}
/// num_rows function for non-MTL types (uses implicit enable_if)
template <typename Collection>
typename traits::num_rows<Collection>::type
inline num_rows(const Collection& c)
{
return traits::num_rows<Collection>()(c);
}
/// num_cols function for non-MTL types (uses implicit enable_if)
template <typename Collection>
typename traits::num_cols<Collection>::type
inline num_cols(const Collection& c)
{
return traits::num_cols<Collection>()(c);
}
/// resize function for vectors
template <typename Collection>
typename traits::num_cols<Collection>::vector_void_type
inline resize(const Collection& c, size_t rows)
{
traits::resize<Collection>()(c, rows);
}
// B(phi)=phi^2*(1-phi)^2 /// resize function for matrices
template <typename Collection>
typename traits::num_cols<Collection>::matrix_void_type
inline resize(const Collection& c, size_t rows, size_t cols)
{
traits::resize<Collection>()(c, rows, cols);
}
template<typename T> template<typename T>
T norm(mtl::dense_vector<T> &b) { return two_norm(b); }; T norm(mtl::dense_vector<T> &b) { return two_norm(b); };
...@@ -185,16 +296,18 @@ namespace vector_operations { ...@@ -185,16 +296,18 @@ namespace vector_operations {
}; };
template<typename Matrix> template<typename Matrix>
double trace(Matrix &M) inline typename boost::enable_if< is_matrix<Matrix>, double >::type
trace(Matrix &M)
{ {
TEST_EXIT(num_rows(M)==num_cols(M))("Matrix dimensions must be equal!\n"); TEST_EXIT(num_rows(M) == num_cols(M))
("Matrix dimensions must be equal!\n");
typename Matrix::value_type sum; typename Matrix::value_type sum;
sum = 0.0; nullify(sum);
for (size_t i = 0; i < num_rows(M); ++i) { for (size_t i = 0; i < num_rows(M); ++i) {
sum += M[i][i]; sum += M[i][i];
} }
return static_cast<double>(sum); return static_cast<double>(sum);
}; }
template<typename T> template<typename T>
Vector<T> mult(const Vector<T>& v1, const Vector<T>& v2) Vector<T> mult(const Vector<T>& v1, const Vector<T>& v2)
...@@ -247,180 +360,137 @@ namespace vector_operations { ...@@ -247,180 +360,137 @@ namespace vector_operations {
return result; return result;
}; };
template<typename T> // copy v -> w
void copy(std::vector<T> &v, WorldVector<T> &w) template<typename Vector1, typename Vector2>
inline void copy(Vector1 &v, Vector2 &w)
{ {
for(int i=0; i<std::min((int)v.size(),(int)w.getSize()); i++) for (size_t i = 0; i < std::min(num_rows(v), num_rows(w)); i++)
w[i]=v[i]; w[i] = v[i];
}; }
template<typename T> // v3:=[v1, v2]
void copy(mtl::dense_vector<T> &v, WorldVector<T> &w) template<typename Vector1, typename Vector2, typename VectorOut>
inline void merge(Vector1 &v1, Vector2 &v2, VectorOut &v3)
{ {
for(int i=0; i<std::min((int)num_rows(v),(int)w.getSize()); i++) resize(v3, num_rows(v1) + num_rows(v2));
w[i]=v[i]; for (size_t i = 0; i < num_rows(v1); i++)
}; v3[i] = v1[i];
for (size_t j = 0; j < num_rows(v2); j++)
template<typename T> v3[j + num_rows(v1)] = v2[j];
void copy(WorldVector<T> &w, mtl::dense_vector<T> &v) }
{
for(int i=0; i<std::min((int)num_rows(v),(int)w.getSize()); i++)
v[i]=w[i];
};
template<typename T> template<typename Vector>
void merge(std::vector<T> &v1, std::vector<T> &v2, std::vector<T> &v3) inline void getMin(const Vector &v, typename ValueType<Vector>::type &minVal, size_t &minIdx)
{ {
v3.resize(v1.size()+v2.size()); typedef typename ValueType<Vector>::type T;
for(int i=0;i<v1.size();i++) TEST_EXIT(num_rows(v) > 0)("getMin of empty vector!\n");
v3[i]=v1[i];
for(int j=0;j<v2.size();j++) minVal = v[0];
v3[j+v1.size()]=v2[j]; minIdx = 0;
}; for (size_t i = 1; i < num_rows(v); i++) {
if (v[i] < minVal) {
template<typename T> minVal = v[i];
void getMin(mtl::dense_vector<T> &v, T &minVal, unsigned &minIdx) minIdx = i;
{
if(num_rows(v)==0) {
minVal = numeric_limits<T>::infinity();
} else {
minVal=v[0]; minIdx=0;
for(unsigned i=1; i<num_rows(v); i++) {
if(v[i]<minVal) { minVal=v[i]; minIdx=i; }
} }
} }
}; }
template<typename T> template<typename Vector>
void getMin(std::vector<T> &v, T &minVal, size_t &minIdx) inline void getMax(const Vector &v, typename ValueType<Vector>::type &maxVal, size_t &maxIdx)
{ {
if (v.size() == 0) { typedef typename ValueType<Vector>::type T;
minVal = numeric_limits<T>::infinity(); TEST_EXIT(num_rows(v) > 0)("getMax of empty vector!\n");
} else {
minVal = v[0]; maxVal = v[0];
minIdx = 0; maxIdx = 0;
for(size_t i = 1; i < v.size(); i++) { for (size_t i = 1; i < num_rows(v); i++) {
if (v[i] < minVal) { if (v[i] > maxVal) {
minVal=v[i]; maxVal = v[i];
minIdx=i; maxIdx = i;
}
} }
} }
}; }
template<typename T> template<typename Vector>
void getMax(mtl::dense_vector<T> &v, T &maxVal, unsigned &maxIdx) inline typename ValueType<Vector>::type
min(const Vector &vec)