diff --git a/AMDiS/CMakeLists.txt b/AMDiS/CMakeLists.txt index 30c4e218089d166b1d4be55345a9f3a952041143..052ed3ea27b9a61b982acf83eb42aa485f1dda34 100644 --- a/AMDiS/CMakeLists.txt +++ b/AMDiS/CMakeLists.txt @@ -685,6 +685,11 @@ INSTALL(FILES ${HEADERS} DESTINATION include/amdis/traits/) list(APPEND deb_add_dirs "include/amdis/traits") +FILE(GLOB HEADERS "${SOURCE_DIR}/utility/*.h*") +INSTALL(FILES ${HEADERS} + DESTINATION include/amdis/utility/) +list(APPEND deb_add_dirs "include/amdis/utility") + FILE(GLOB HEADERS "${SOURCE_DIR}/time/*.h") INSTALL(FILES ${HEADERS} DESTINATION include/amdis/time/) diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc index 1994118b9c53354df89a6698acac33950b1c3fee..687b8e88e99731452f628eb612917d960df15699 100644 --- a/AMDiS/src/DOFVector.cc +++ b/AMDiS/src/DOFVector.cc @@ -486,7 +486,7 @@ namespace AMDiS { elInfo2->worldToCoord(worldVec, &coords2); bool isPositive = true; - for (int j = 0; j < coords2.size(); j++) { + for (int j = 0; j < coords2.getSize(); j++) { if (coords2[j] < -0.00001) { isPositive = false; break; diff --git a/AMDiS/src/ElInfo.h b/AMDiS/src/ElInfo.h index f4dc850f70bf98cbc18c282ed31b1289c89f2479..14ce45288352f12d06b646d7441710b1a07cebc6 100644 --- a/AMDiS/src/ElInfo.h +++ b/AMDiS/src/ElInfo.h @@ -215,8 +215,8 @@ namespace AMDiS { inline void getGrdLambda(mtl::dense2D<double>& grd_lam) { - grd_lam.change_dim(grdLambda.size(), Global::getGeo(WORLD)); - for (size_t i = 0; i < static_cast<size_t>(grdLambda.size()); i++) + grd_lam.change_dim(grdLambda.getSize(), Global::getGeo(WORLD)); + for (size_t i = 0; i < static_cast<size_t>(grdLambda.getSize()); i++) for (size_t j = 0; j < static_cast<size_t>(Global::getGeo(WORLD)); j++) grd_lam(i,j) = grd_lam[i][j]; } diff --git a/AMDiS/src/FixVec.h b/AMDiS/src/FixVec.h index 8a06436d624ade707bf7dbfed4a8dc103bec839e..b5a3346bf8132f8fe3377acf7fa073c4a97b2c7f 100644 --- a/AMDiS/src/FixVec.h +++ b/AMDiS/src/FixVec.h @@ -53,6 +53,8 @@ namespace AMDiS { template<typename T,GeoIndex d> class FixVec : public Vector<T> { + typedef FixVec self; + typedef Vector<T> super; public: /// Constructor without initialisation. initType must be NO_INIT. If dim is @@ -80,6 +82,10 @@ namespace AMDiS { TEST_EXIT_DBG(initType == DEFAULT_VALUE)("wrong initType or wrong initializer\n"); this->set(ini); } + + FixVec(self const& other) + : super(other) + { } /// Initialisation for dim. inline void init(int dim) @@ -94,10 +100,10 @@ namespace AMDiS { } /// Returns the \ref size_ of the FixVec. - inline int size() const - { - return this->getSize(); - } +// inline int size() const +// { +// return this->getSize(); +// } protected: /// Determines needed vector size. @@ -178,7 +184,7 @@ namespace AMDiS { } /// Destructor - virtual ~VectorOfFixVecs() + ~VectorOfFixVecs() { for (int i = 0; i < size; i++) delete vec[i]; @@ -274,7 +280,7 @@ namespace AMDiS { } /// destructor - virtual ~MatrixOfFixVecs() + ~MatrixOfFixVecs() { for (VectorOfFixVecs<FixVecType>** i = &vec[0]; i < &vec[rows]; i++) delete *i; @@ -417,18 +423,19 @@ namespace AMDiS { {} /// Copy assignement operator - self& operator=(self other) - { - swap(*this, other); - return *this; - } +// self& operator=(self const& other) +// { +// assert( Global::getGeo(WORLD) == other.getSize() ); +// this->setValues(other.getValArray()); +// return *this; +// } /// Assignement operator template <typename S> self& operator=(const Vector<S>& other) { - TEST_EXIT_DBG( other.getSize() == Global::getGeo(WORLD) ) + TEST_EXIT_DBG( other.getSize() == this->size ) ("Wrong dimensions in assignment.\n"); this->setValues(other.getValArray()); return *this; @@ -503,26 +510,27 @@ namespace AMDiS { { } /// Copy assignment operator - self& operator=(self other) - { - swap(*this, other); - return *this; - } +// self& operator=(self other) +// { +// swap(*this, other); +// return *this; +// } /// Assignment operator template <typename S> self& operator=(const Matrix<S>& other) { - TEST_EXIT_DBG( other.getNumOfRows() == Global::getGeo(WORLD) && - other.getNumOfCols() == Global::getGeo(WORLD) ) - ("Wrong dimensions in assignment.\n"); + TEST_EXIT_DBG( this->size == other.getSize() && + this->rows == other.getNumRows() && + this->cols == other.getNumCols() ) + ("Wrong dimensions in assignment.\n"); this->setValues(other.getValArray()); return *this; } /// Assignment operator for scalars template <typename S> - typename enable_if<boost::is_convertible<S, T>, WorldMatrix<T> >::type & + typename enable_if<boost::is_convertible<S, T>, self>::type & operator=(S value) { this->set(value); diff --git a/AMDiS/src/FixVec.hh b/AMDiS/src/FixVec.hh index fcad92748e653269798f323b1122584495129d5e..d644c64b9d2b821744d72f3c354731f4062e9f30 100644 --- a/AMDiS/src/FixVec.hh +++ b/AMDiS/src/FixVec.hh @@ -87,8 +87,8 @@ namespace AMDiS { { FUNCNAME_DBG("WorldMatrix<T>::vecProduct()"); - TEST_EXIT_DBG(v1.getSize() == v2.getSize())("invalid size 1\n"); - TEST_EXIT_DBG(v1.getSize() == this->getSize())("invalid size 2\n"); + TEST_EXIT_DBG(v1.getSize() == v2.getSize())("size(v1) != size(v2), %d != %d\n", v1.getSize(), v2.getSize()); + TEST_EXIT_DBG(v1.getSize() == this->getNumRows())("size(v1) != num_rows(this), %d != %d\n", v1.getSize(), this->getNumRows()); T *thisIt = this->begin(); diff --git a/AMDiS/src/Global.cc b/AMDiS/src/Global.cc index 3457a8a0d29c8b4594c2a257aa11fce5dd22a2ca..4b433c0d54c5b6d6ff712db877c1418c2844d3fd 100644 --- a/AMDiS/src/Global.cc +++ b/AMDiS/src/Global.cc @@ -61,7 +61,7 @@ namespace AMDiS { FUNCNAME("Msg::wait()"); if (w) { - std::string line; + char line; MSG("wait for <enter> ..."); std::cin >> line; // char* result = fgets(line, 9, stdin); diff --git a/AMDiS/src/Global.h b/AMDiS/src/Global.h index 2f44fdc50c4b19a315bf2039f165afa83fdbed67..d49cc3b2352de9b92d3d40c8c9d1ce03fa4b5193 100644 --- a/AMDiS/src/Global.h +++ b/AMDiS/src/Global.h @@ -361,7 +361,7 @@ namespace AMDiS { /// function name nn for message output via MSG, WARNING, ... #define FUNCNAME(nn) const char *funcName; funcName = nn; -#if (DEBUG == 0) +#ifdef NDEBUG #define FUNCNAME_DBG(nn) #else #define FUNCNAME_DBG(nn) const char *funcName; funcName = nn; @@ -386,18 +386,18 @@ namespace AMDiS { #define TEST_EXIT(test) if ((test));else ERROR_EXIT /// In debug mode, it corresponds to ERROR_EXIT, otherwise it is noop. -#if (DEBUG == 0) +#ifdef NDEBUG #define TEST_EXIT_DBG(test) if (false) Msg::catch_error_exit #define DBG_VAR(var) #else -#define TEST_EXIT_DBG(test) if ((test));else ERROR_EXIT + #define TEST_EXIT_DBG(test) if ((test));else ERROR_EXIT #define DBG_VAR(var) var #endif /// prints a message #define MSG Msg::print_funcname(funcName), Msg::print -#if (DEBUG == 0) +#ifdef NDEBUG #define MSG_DBG #else #define MSG_DBG Msg::print_funcname(funcName), Msg::print diff --git a/AMDiS/src/MatrixVector.h b/AMDiS/src/MatrixVector.h index 3098d9f5eb1da62a00e517f14333c987f4e0d4d2..96fefd1a98b16161a65890d7a60e687c7a108d61 100644 --- a/AMDiS/src/MatrixVector.h +++ b/AMDiS/src/MatrixVector.h @@ -56,31 +56,11 @@ namespace AMDiS { { setValues(other.valArray); } - - /// Copy constructor for other of different value_type - // TODO: notwendig? -// template <typename S> -// Vector(Vector<S> const& rhs) -// : Serializable(), -// size(rhs.getSize()), -// valArray(size ? new T[size] : NULL) -// { -// setValues(rhs.valArray); -// } - -// #if HAS_RVALUE_REFERENCES -// /// move constructor -// Vector(self&& other) -// : Vector() -// { -// swap(*this, other); -// } -// #endif /// Destructor. - virtual ~Vector() + ~Vector() { - if (valArray != NULL) { + if (valArray) { delete [] valArray; valArray = NULL; } @@ -95,7 +75,7 @@ namespace AMDiS { inline void resize(int s) { if (size != s && s > 0) { - if (valArray != NULL) + if (valArray) delete [] valArray; valArray = new T[s]; size = s; @@ -103,24 +83,9 @@ namespace AMDiS { } /// Assignement operator - template <typename S> - typename disable_if< boost::is_same<S, T>, self >::type & - operator=(Vector<S> const& other) - { - assert( size == other.getSize() ); -// resize(other.getSize()); - setValues(other.getValArray()); - return *this; - } - - /// copy assignment operator -// self& operator=(self other) -// { -// swap(*this, other); - self& operator=(self const& other) { - assert( size == other.getSize() ); + resize( other.getSize() ); setValues(other.getValArray()); return *this; } @@ -140,73 +105,33 @@ namespace AMDiS { setValues(vec); return *this; } - - /// A swap function for Vector, used in the Copy-and-swap idiom - // need non-templated arguments in order to eliminate a friend declaration warning in gcc - friend void swap(Vector& first, Vector& second) - { - using std::swap; // enable ADL - swap(first.size, second.size); - swap(first.valArray, second.valArray); - } /// Sets all entries to scal. template <typename S> typename enable_if< boost::is_convertible<S, T> >::type set(S value) { - for (T *thisIt = this->begin(); thisIt != this->end(); ++thisIt) - *thisIt = value; -// std::fill(begin(), end(), value); + std::fill(begin(), end(), value); } /// Sets all entries. template <typename S> inline void setValues(const S* values) { - T *thisIt; - const S *valuesIt; - for (thisIt = this->begin(), valuesIt = values; - thisIt != this->end(); - ++thisIt, ++valuesIt) - *thisIt = *valuesIt; -// std::copy(values, values + size, valArray); - } - - /// Comparison operator. - inline bool operator==(Vector<T> const& rhs) const - { - if (size != rhs.size) - return false; - - T const* rhsIt; - T const* thisIt; - for (rhsIt = rhs.begin(), thisIt = this->begin(); - rhsIt != rhs.end(); - ++rhsIt, ++thisIt) - if (*thisIt != *rhsIt) - return false; - - return true; - } - - /// Comparison operator. - inline bool operator!=(Vector<T> const& rhs) const - { - return !(*this==rhs); + std::copy(values, values + size, begin()); } /// Access to the i-th vector element. inline T& operator[](int i) { -// TEST_EXIT_DBG(i < size && i >= 0)("Invalid index %d!\n", i); + TEST_EXIT_DBG((unsigned)i < (unsigned)size)("Invalid index %d!\n", i); return valArray[i]; } /// Access to the i-th vector element for const vectors. inline const T& operator[] (int i) const { -// TEST_EXIT_DBG(i < size && i >= 0)("Invalid index %d!\n", i); + TEST_EXIT_DBG((unsigned)i < (unsigned)size)("Invalid index %d!\n", i); return valArray[i]; } @@ -288,15 +213,6 @@ namespace AMDiS { rows(other.getNumRows()), cols(other.getNumCols()) { } - -// #if HAS_RVALUE_REFERENCES -// /// move constructor -// Matrix(self&& other) -// : Matrix() -// { -// swap(*this, other); -// } -// #endif /// Changes the size of the matrix to newRows x newCols. inline void resize(int newRows, int newCols) @@ -307,60 +223,22 @@ namespace AMDiS { cols = newCols; } } - -// self& operator=(self other) -// { -// swap(*this, other); -// return *this; -// } - + self& operator=(self const& other) { - assert( this->size == other.getSize() && rows == other.getNumRows() && cols == other.getNumCols() ); - this->setValues(other.getValArray()); - return *this; - } - - /// Assignement operator - template <typename S> - inline self& operator=(const Matrix<S>& other) - { - assert( rows == other.getNumRows() && cols == other.getNumCols() ); -// resize(other.getNumRows(), other.getNumCols()); + resize( other.getNumRows(), other.getNumCols() ); this->setValues(other.getValArray()); return *this; } /// Assignement operator for scalars. - self& operator=(T value) + template <typename S> + typename enable_if< boost::is_convertible<S, T>, self>::type & + operator=(S value) { - super::set(value); + this->set(value); return *this; } - - /// A swap function for Matrix, used in the Copy-and-swap idiom - // need non-templated arguments in order to eliminate a friend declaration warning in gcc - friend void swap(Matrix& first, Matrix& second) - { - using std::swap; // enable ADL - swap(first.rows, second.rows); - swap(first.cols, second.cols); - swap(static_cast<super&>(first), static_cast<super&>(second)); - } - - /// Comparison operator. - inline bool operator==(const self& rhs) const - { - if (rows != rhs.getNumRows()) return false; - if (cols != rhs.getNumCols()) return false; - return super::operator==(rhs); - } - - /// Comparison operator. - inline bool operator!=(const self& rhs) const - { - return !(*this == rhs); - } /// Access to i-th matrix row. inline T *operator[](int i) diff --git a/AMDiS/src/MatrixVectorOperations.h b/AMDiS/src/MatrixVectorOperations.h index f3bfb1d9913627f975be28f533e79623cfbf5861..45cfeca74f8028320cbed59908a374ac1c71d4dc 100644 --- a/AMDiS/src/MatrixVectorOperations.h +++ b/AMDiS/src/MatrixVectorOperations.h @@ -473,6 +473,20 @@ namespace AMDiS { return m1; } + /// matrix += matrix + template <typename T, typename S> + Matrix<T>& operator+=(Matrix<T>& m1, const Matrix<S>& m2) + { + T* m1It; + S const* m2It; + for (m1It = m1.begin(), m2It = m2.begin(); + m1It != m1.end(); + m1It++, m2It++) + *m1It += *m2It; + + return m1; + } + /// matrix := matrix + matrix template <typename T, typename S> WorldMatrix<T> operator+(WorldMatrix<T> M1, const WorldMatrix<S>& M2 ) @@ -480,6 +494,14 @@ namespace AMDiS { M1 += M2; return M1; } + + /// matrix := matrix + matrix + template <typename T, typename S> + Matrix<T> operator+(Matrix<T> M1, const Matrix<S>& M2 ) + { + M1 += M2; + return M1; + } /// matrix -= matrix template <typename T, typename S> @@ -494,6 +516,20 @@ namespace AMDiS { return m1; } + + /// matrix -= matrix + template <typename T, typename S> + Matrix<T>& operator-=(Matrix<T>& m1, const Matrix<S>& m2) + { + T *m1It; + S const* m2It; + for (m1It = m1.begin(), m2It = m2.begin(); + m1It != m1.end(); + m1It++, m2It++) + *m1It -= *m2It; + + return m1; + } /// matrix := matrix - matrix template <typename T, typename S> @@ -502,6 +538,14 @@ namespace AMDiS { M1 -= M2; return M1; } + + /// matrix := matrix - matrix + template <typename T, typename S> + Matrix<T> operator-(Matrix<T> M1, const Matrix<S>& M2 ) + { + M1 -= M2; + return M1; + } // unary minus operators // --------------------- @@ -548,6 +592,50 @@ namespace AMDiS { return true; } + + /// Comparison operator for Vector<T>. + template <typename T> + inline bool operator==(Vector<T> const& lhs, Vector<T> const& rhs) + { + if (lhs.getSize() != rhs.getSize()) + return false; + + T const* rhsIt; + T const* lhsIt; + for (rhsIt = rhs.begin(), lhsIt = lhs.begin(); + rhsIt != rhs.end(); + ++rhsIt, ++lhsIt) + if (*lhsIt != *rhsIt) + return false; + + return true; + } + + /// Comparison operator for Vector<T>. + template <typename T> + inline bool operator!=(Vector<T> const& lhs, Vector<T> const& rhs) + { + return !(lhs == rhs); + } + + + + /// Comparison operator for Matrix<T>. + template <typename T> + inline bool operator==(Matrix<T> const& lhs, Matrix<T> const& rhs) + { + if (lhs.getNumRows() != rhs.getNumRows()) return false; + if (lhs.getNumCols() != rhs.getNumCols()) return false; + return (static_cast<Vector<T> const&>(lhs) == static_cast<Vector<T> const&>(rhs)); + } + + /// Comparison operator for Matrix<T>. + template <typename T> + inline bool operator!=(Matrix<T> const& lhs, Matrix<T> const& rhs) + { + return !(lhs == rhs); + } + // special operators // ----------------- diff --git a/AMDiS/src/Recovery.h b/AMDiS/src/Recovery.h index 468e601074e4fb3fe691526107b282f35ec545e1..cf637ac239bb5b59b1641c14d267d91c22d4b0ff 100644 --- a/AMDiS/src/Recovery.h +++ b/AMDiS/src/Recovery.h @@ -42,7 +42,7 @@ namespace AMDiS { exponent(expon) { degree_ = exponent[0]; - for (int i = 1; i<exponent.size(); i++) + for (int i = 1; i<exponent.getSize(); i++) degree_ += exponent[i]; } @@ -52,7 +52,7 @@ namespace AMDiS { const WorldVector<double>& z) const { double result = std::pow(y[0] - z[0], double(exponent[0])); - for (int i = 1; i < exponent.size(); i++) + for (int i = 1; i < exponent.getSize(); i++) result *= std::pow(y[i] - z[i], double(exponent[i])); return result; diff --git a/AMDiS/src/VertexInfo.h b/AMDiS/src/VertexInfo.h index 7d6756b885f20aad63c178fb52e1ddb9937074fb..84db86d56209a8727b96a14e12f3e92d84d85657 100644 --- a/AMDiS/src/VertexInfo.h +++ b/AMDiS/src/VertexInfo.h @@ -26,6 +26,7 @@ #define AMDIS_VERTEXINFO_H #include "FixVec.h" +#include "MatrixVectorOperations.h" namespace AMDiS { diff --git a/AMDiS/src/expressions/value_expr.hpp b/AMDiS/src/expressions/value_expr.hpp index aa9277b33eb526f14cf473ab7a875b03c0e97b64..8bff587e1931c1dc51b0c3a973b84cc6ce226f46 100644 --- a/AMDiS/src/expressions/value_expr.hpp +++ b/AMDiS/src/expressions/value_expr.hpp @@ -164,11 +164,10 @@ namespace AMDiS struct Eye : public LazyOperatorTermBase { typedef Matrix<int> value_type; - value_type M; size_t n; + value_type M; - Eye(int n) : M(n,n) { - assert((N >= 0 && n < 0) || (N < 0 && n >= 0)); + Eye(int n_) : n(n_ < 0 ? N : n_), M(n,n) { M = 0; for (size_t i = 0; i < n; ++i) M[i][i] = 1; diff --git a/AMDiS/src/io/VtkVectorWriter.cc b/AMDiS/src/io/VtkVectorWriter.cc index 817520876970f945d44381915b33d9605e82703e..baa31816a514460814bfb4ea8cd00cb0968ba0f2 100644 --- a/AMDiS/src/io/VtkVectorWriter.cc +++ b/AMDiS/src/io/VtkVectorWriter.cc @@ -84,18 +84,18 @@ namespace AMDiS { DOFVector<WorldVector<double> > *newValues = new DOFVector<WorldVector<double> >(values[0]->getFeSpace(), values[0]->getName()); WorldVector<DOFIterator<double>* > iterators; - for (size_t i = 0; i < static_cast<size_t>(values.getSize()); i++) + for (int i = 0; i < values.getSize(); i++) iterators[i] = new DOFIterator<double>(values[i],USED_DOFS); - for (size_t i = 0; i < static_cast<size_t>(iterators.getSize()); i++) + for (int i = 0; i < iterators.getSize(); i++) iterators[i]->reset(); DOFIterator<WorldVector<double> > resultIter(newValues, USED_DOFS); for(resultIter.reset(); !resultIter.end(); resultIter++) { - for (size_t i = 0; i < static_cast<size_t>(iterators.getSize()); i++) + for (int i = 0; i < iterators.getSize(); i++) (*resultIter)[i] = *(*(iterators[i])); - for (size_t i = 0; i < static_cast<size_t>(iterators.size()); i++) + for (int i = 0; i < iterators.getSize(); i++) (*(iterators[i]))++; } diff --git a/AMDiS/src/io/detail/Arh2Reader.cc b/AMDiS/src/io/detail/Arh2Reader.cc index 7a3e620f61a423b90d655213f64312726f41faeb..7e12a074db5f68ce439fbdb9c88c8e2084320e2b 100644 --- a/AMDiS/src/io/detail/Arh2Reader.cc +++ b/AMDiS/src/io/detail/Arh2Reader.cc @@ -323,14 +323,14 @@ namespace AMDiS { namespace io { continue; } int j; - for(j = 0; j < nDOF->size(); j++) + for(j = 0; j < nDOF->getSize(); j++) { if((*nDOF)[j] != feSpaceDOFs[vecsFeSpaceNum[i]][j]) { break; } } - if(j == nDOF->size()) + if(j == nDOF->getSize()) { vecs[i] = tmpVecs[k]; break; @@ -351,7 +351,7 @@ namespace AMDiS { namespace io { ("Vecs number %i has no DOFAdmin. Should not happen.\n", i); DimVec<int>* nDOF = vecs[i]->getFeSpace()->getBasisFcts()->getNumberOfDofs(); - for(int j = 0; j < nDOF->size(); j++) + for(int j = 0; j < nDOF->getSize(); j++) { TEST_EXIT((*nDOF)[j] == feSpaceDOFs[vecsFeSpaceNum[i]][j]) ("The fespace of vec number %i is not equal to the correspond fespace.\n", i+1); diff --git a/AMDiS/src/io/detail/Arh2Writer.cc b/AMDiS/src/io/detail/Arh2Writer.cc index 5b79c0162715bc508afcbcf666acfe973973203a..6e95ce8e4248cf2cdbeccf4ab7ce397d7813807c 100644 --- a/AMDiS/src/io/detail/Arh2Writer.cc +++ b/AMDiS/src/io/detail/Arh2Writer.cc @@ -188,12 +188,12 @@ namespace AMDiS { namespace io { for(feSpaceIt = feSpaces.begin(); feSpaceIt != feSpaces.end(); feSpaceIt++, i++) { DimVec<int>* nDOF = feSpaceIt->first->getBasisFcts()->getNumberOfDofs(); - for(int j = 0; j < nDOF->size(); j++) + for(int j = 0; j < nDOF->getSize(); j++) { posDOFs = (*nDOF)[j]; file.write(reinterpret_cast<char*>(&posDOFs), 4); } - for(size_t j = nDOF->size(); j < 4 ; j++) + for(size_t j = nDOF->getSize(); j < 4 ; j++) { posDOFs = 0; file.write(reinterpret_cast<char*>(&posDOFs), 4); diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc index 0d5909a2082344cad7461cbe80d1185d0dcca249..dbb2bfa527cfab86d9354c9b4a0a08d90bb3be4e 100644 --- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc +++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc @@ -159,7 +159,7 @@ namespace AMDiS { namespace Parallel { initSolver(kspInterior); KSPGetPC(kspInterior, &pcInterior); - initPreconditioner(pcInterior); + initPreconditioner(*seqMat, mat[0][0]); #if (DEBUG != 0) diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.h b/AMDiS/src/parallel/PetscSolverGlobalMatrix.h index a2a9a0343fb02da1d243fbbaeea02f776a0ce168..aee818fc56adcace4f81976ece3b7fd8be3392b0 100644 --- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.h +++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.h @@ -113,6 +113,10 @@ namespace AMDiS virtual void exitSolver(KSP &ksp); virtual void initPreconditioner(PC pc); + virtual void initPreconditioner(const Matrix<DOFMatrix*>& A, const Mat& fullMatrix) + { + initPreconditioner(getPc()); + } virtual void exitPreconditioner(PC pc); diff --git a/AMDiS/src/parallel/StdMpi.cc b/AMDiS/src/parallel/StdMpi.cc index 15ae83026f73d4f78d76477fb09c670697383125..ba7ac5d9589c145037638f489bc3cb7b7f5a4526 100644 --- a/AMDiS/src/parallel/StdMpi.cc +++ b/AMDiS/src/parallel/StdMpi.cc @@ -37,7 +37,12 @@ namespace AMDiS { namespace Parallel { MPI_Datatype StdMpiHelper<map<BoundaryType, map<DegreeOfFreedom, DegreeOfFreedom> > >::mpiDataType = MPI_INT; MPI_Datatype StdMpiHelper<vector<std::pair<int, int> > >::mpiDataType = MPI_INT; MPI_Datatype StdMpiHelper<vector<WorldVector<double> > >::mpiDataType = MPI_DOUBLE; + MPI_Datatype StdMpiHelper<vector<WorldVector<int> > >::mpiDataType = MPI_INT; + MPI_Datatype StdMpiHelper<vector<Vector<int> > >::mpiDataType = MPI_INT; MPI_Datatype StdMpiHelper<vector<WorldMatrix<double> > >::mpiDataType = MPI_DOUBLE; + MPI_Datatype StdMpiHelper<vector<WorldVector<WorldVector<double> > > >::mpiDataType = MPI_DOUBLE; + MPI_Datatype StdMpiHelper<vector<WorldMatrix<int> > >::mpiDataType = MPI_INT; + MPI_Datatype StdMpiHelper<vector<Matrix<int> > >::mpiDataType = MPI_INT; MPI_Datatype StdMpiHelper<map<WorldVector<double>, int> >::mpiDataType = MPI_DOUBLE; MPI_Datatype StdMpiHelper<map<int, WorldVector<double> > >::mpiDataType = MPI_DOUBLE; MPI_Datatype StdMpiHelper<vector<vector<WorldVector<double> > > >::mpiDataType = MPI_DOUBLE; @@ -468,6 +473,63 @@ namespace AMDiS { namespace Parallel { + // T = vector<WorldVector<int> > + + int StdMpiHelper<vector<WorldVector<int> > >::getBufferSize(vector<WorldVector<int> > &data) + { + return data.size() * Global::getGeo(WORLD); + } + + void StdMpiHelper<vector<WorldVector<int> > >::createBuffer(vector<WorldVector<int> > &data, int *buf) + { + int dimOfWorld = Global::getGeo(WORLD); + int pos = 0; + for (unsigned int i = 0; i < data.size(); i++) + for (int j = 0; j < dimOfWorld; j++) + buf[pos++] = data[i][j]; + } + + void StdMpiHelper<vector<WorldVector<int> > >::makeFromBuffer(vector<WorldVector<int> > &data, int *buf, int bufSize) + { + int dimOfWorld = Global::getGeo(WORLD); + TEST_EXIT(bufSize % Global::getGeo(WORLD) == 0)("This should not happen!\n"); + + int pos = 0; + data.resize(bufSize / Global::getGeo(WORLD)); + for (unsigned int i = 0; i < data.size(); i++) + for (int j = 0; j < dimOfWorld; j++) + data[i][j] = buf[pos++]; + } + + + + // T = vector<Vector<int> > + + int StdMpiHelper<vector<Vector<int> > >::getBufferSize(vector<Vector<int> > &data) + { + return data.size() > 0 ? data.size() * data[0].getSize() : 0; + } + + void StdMpiHelper<vector<Vector<int> > >::createBuffer(vector<Vector<int> > &data, int *buf) + { + int pos = 0; + buf[pos++] = data.size() > 0 ? data[0].getSize() : 0; + for (unsigned int i = 0; i < data.size(); i++) + for (int j = 0; j < data[i].getSize(); j++) + buf[pos++] = data[i][j]; + } + + void StdMpiHelper<vector<Vector<int> > >::makeFromBuffer(vector<Vector<int> > &data, int *buf, int bufSize) + { + int pos = 0; + int size = buf[pos++]; + data.resize(bufSize / size); + for (unsigned int i = 0; i < data.size(); i++) + for (int j = 0; j < size; j++) + data[i][j] = buf[pos++]; + } + + // T = vector<WorldMatrix<double> > int StdMpiHelper<vector<WorldMatrix<double> > >::getBufferSize(vector<WorldMatrix<double> > &data) @@ -499,6 +561,103 @@ namespace AMDiS { namespace Parallel { } + // T = vector<WorldMatrix<double> > + + int StdMpiHelper<vector<WorldVector<WorldVector<double> > > >::getBufferSize(vector<WorldVector<WorldVector<double> > > &data) + { + return data.size() * sqr(Global::getGeo(WORLD)); + } + + void StdMpiHelper<vector<WorldVector<WorldVector<double> > > >::createBuffer(vector<WorldVector<WorldVector<double> > > &data, double *buf) + { + int dimOfWorld = Global::getGeo(WORLD); + int pos = 0; + for (unsigned int i = 0; i < data.size(); i++) + for (int j = 0; j < dimOfWorld; j++) + for (int k = 0; k < dimOfWorld; k++) + buf[pos++] = data[i][j][k]; + } + + void StdMpiHelper<vector<WorldVector<WorldVector<double> > > >::makeFromBuffer(vector<WorldVector<WorldVector<double> > > &data, double *buf, int bufSize) + { + int dimOfWorld = Global::getGeo(WORLD); + TEST_EXIT(bufSize % sqr(Global::getGeo(WORLD)) == 0)("This should not happen!\n"); + + int pos = 0; + data.resize(bufSize / sqr(dimOfWorld)); + for (unsigned int i = 0; i < data.size(); i++) + for (int j = 0; j < dimOfWorld; j++) + for (int k = 0; k < dimOfWorld; k++) + data[i][j][k] = buf[pos++]; + } + + + + // T = vector<WorldMatrix<int> > + + int StdMpiHelper<vector<WorldMatrix<int> > >::getBufferSize(vector<WorldMatrix<int> > &data) + { + return data.size() * sqr(Global::getGeo(WORLD)); + } + + void StdMpiHelper<vector<WorldMatrix<int> > >::createBuffer(vector<WorldMatrix<int> > &data, int *buf) + { + int dimOfWorld = Global::getGeo(WORLD); + int pos = 0; + for (unsigned int i = 0; i < data.size(); i++) + for (int j = 0; j < dimOfWorld; j++) + for (int k = 0; k < dimOfWorld; k++) + buf[pos++] = data[i][j][k]; + } + + void StdMpiHelper<vector<WorldMatrix<int> > >::makeFromBuffer(vector<WorldMatrix<int> > &data, int *buf, int bufSize) + { + int dimOfWorld = Global::getGeo(WORLD); + TEST_EXIT(bufSize % sqr(Global::getGeo(WORLD)) == 0)("This should not happen!\n"); + + int pos = 0; + data.resize(bufSize / sqr(dimOfWorld)); + for (unsigned int i = 0; i < data.size(); i++) + for (int j = 0; j < dimOfWorld; j++) + for (int k = 0; k < dimOfWorld; k++) + data[i][j][k] = buf[pos++]; + } + + + // T = vector<Matrix<int> > + + int StdMpiHelper<vector<Matrix<int> > >::getBufferSize(vector<Matrix<int> > &data) + { + return data.size() > 0 ? data.size() * data[0].getNumRows() * data[0].getNumCols() + 2 : 0; + } + + void StdMpiHelper<vector<Matrix<int> > >::createBuffer(vector<Matrix<int> > &data, int *buf) + { + int dimOfWorld = Global::getGeo(WORLD); + int pos = 0; + buf[pos++] = data.size() > 0 ? data[0].getNumRows() : 0; + buf[pos++] = data.size() > 0 ? data[0].getNumCols() : 0; + for (unsigned int i = 0; i < data.size(); i++) + for (int j = 0; j < data[i].getNumRows(); j++) + for (int k = 0; k < data[i].getNumCols(); k++) + buf[pos++] = data[i][j][k]; + } + + void StdMpiHelper<vector<Matrix<int> > >::makeFromBuffer(vector<Matrix<int> > &data, int *buf, int bufSize) + { + int pos = 0; + int n_rows = buf[pos++]; + int n_cols = buf[pos++]; + data.resize(bufSize / (n_rows * n_cols)); + for (unsigned int i = 0; i < data.size(); i++) { + data[i].resize(n_rows, n_cols); + for (int j = 0; j < n_rows; j++) + for (int k = 0; k < n_cols; k++) + data[i][j][k] = buf[pos++]; + } + } + + // T = map<WorldVector<double>, int> diff --git a/AMDiS/src/parallel/StdMpi.h b/AMDiS/src/parallel/StdMpi.h index 9fb376566357260591ff35110dcc5215149edf99..a9258df349914330e27417f74a3e7736e7ed96af 100644 --- a/AMDiS/src/parallel/StdMpi.h +++ b/AMDiS/src/parallel/StdMpi.h @@ -226,6 +226,32 @@ namespace AMDiS { namespace Parallel { static void makeFromBuffer(std::vector<WorldVector<double> > &data, double *buf, int bufSize); }; + template<> + struct StdMpiHelper<std::vector<WorldVector<int> > > { + static MPI_Datatype mpiDataType; + + typedef int cppDataType; + + static int getBufferSize(std::vector<WorldVector<int> > &data); + + static void createBuffer(std::vector<WorldVector<int> > &data, int *buf); + + static void makeFromBuffer(std::vector<WorldVector<int> > &data, int *buf, int bufSize); + }; + + template<> + struct StdMpiHelper<std::vector<Vector<int> > > { + static MPI_Datatype mpiDataType; + + typedef int cppDataType; + + static int getBufferSize(std::vector<Vector<int> > &data); + + static void createBuffer(std::vector<Vector<int> > &data, int *buf); + + static void makeFromBuffer(std::vector<Vector<int> > &data, int *buf, int bufSize); + }; + template<> struct StdMpiHelper<std::vector<WorldMatrix<double> > > { static MPI_Datatype mpiDataType; @@ -238,6 +264,45 @@ namespace AMDiS { namespace Parallel { static void makeFromBuffer(std::vector<WorldMatrix<double> > &data, double *buf, int bufSize); }; + + template<> + struct StdMpiHelper<std::vector<WorldMatrix<int> > > { + static MPI_Datatype mpiDataType; + + typedef int cppDataType; + + static int getBufferSize(std::vector<WorldMatrix<int> > &data); + + static void createBuffer(std::vector<WorldMatrix<int> > &data, int *buf); + + static void makeFromBuffer(std::vector<WorldMatrix<int> > &data, int *buf, int bufSize); + }; + + template<> + struct StdMpiHelper<std::vector<WorldVector<WorldVector<double> > > > { + static MPI_Datatype mpiDataType; + + typedef double cppDataType; + + static int getBufferSize(std::vector<WorldVector<WorldVector<double> > > &data); + + static void createBuffer(std::vector<WorldVector<WorldVector<double> > > &data, double *buf); + + static void makeFromBuffer(std::vector<WorldVector<WorldVector<double> > > &data, double *buf, int bufSize); + }; + + template<> + struct StdMpiHelper<std::vector<Matrix<int> > > { + static MPI_Datatype mpiDataType; + + typedef int cppDataType; + + static int getBufferSize(std::vector<Matrix<int> > &data); + + static void createBuffer(std::vector<Matrix<int> > &data, int *buf); + + static void makeFromBuffer(std::vector<Matrix<int> > &data, int *buf, int bufSize); + }; // ParallelDebug::testGlobalIndexByCoords::CoordsIndexMap diff --git a/AMDiS/src/reinit/BoundaryElementDist.cc b/AMDiS/src/reinit/BoundaryElementDist.cc index d3f6c7b205853a7afca5f2f4d5ed344f1c18a7bc..4bf66a353bf9c005538e2d935fccaf29912639c2 100644 --- a/AMDiS/src/reinit/BoundaryElementDist.cc +++ b/AMDiS/src/reinit/BoundaryElementDist.cc @@ -48,7 +48,7 @@ BoundaryElementDist::calcNormal_2d( { FUNCNAME("BoundaryElementDist::calcNormal_2d"); - TEST_EXIT(vecs.size() == dim)("illegal number of world vectors in vecs !\n"); + TEST_EXIT(vecs.getSize() == dim)("illegal number of world vectors in vecs !\n"); double norm2 = 0.0; double val; @@ -71,7 +71,7 @@ BoundaryElementDist::calcNormal_3d( { FUNCNAME("BoundaryElementDist::calcNormal_3d"); - TEST_EXIT(vecs.size() == dim)("illegal number of world vectors in vecs !\n"); + TEST_EXIT(vecs.getSize() == dim)("illegal number of world vectors in vecs !\n"); WorldVector<double> A = vecs[1]-vecs[0]; WorldVector<double> B = vecs[2]-vecs[0]; diff --git a/AMDiS/src/solver/Mapper.h b/AMDiS/src/solver/Mapper.h index 0078f1d2335c8886be68fa14939649a62de4b240..d32f206916a14460185fffc3d93045b50e732359 100644 --- a/AMDiS/src/solver/Mapper.h +++ b/AMDiS/src/solver/Mapper.h @@ -44,6 +44,11 @@ namespace AMDiS { { typedef MTLTypes::size_type size_type; + template <class BaseInserter> + struct Inserter { + typedef mtl::matrix::mapped_inserter<BaseInserter, Derived> type; + }; + /// set the current block row inline void setRow( unsigned int r) { self().setRow(r); } @@ -92,6 +97,37 @@ namespace AMDiS { { typedef unsigned int size_type; + /// Default constructor + BlockMapper() + : nComp(0), rowOffset(0), colOffset(0), nrow(0), ncol(0), sizes(0) + {} + + BlockMapper(BlockMapper const& other) + : nComp(other.nComp), + rowOffset(other.rowOffset), + colOffset(other.colOffset), + nrow(other.nrow), + ncol(other.ncol) + { + sizes.resize(nComp); + std::copy(other.sizes.begin(), other.sizes.end(), sizes.begin()); + } + + BlockMapper& operator=(BlockMapper const& other) + { + nComp = other.getNumComponents(); + rowOffset = other.row(0); + colOffset = other.col(0); + nrow = other.getNumRows(); + ncol = other.getNumCols(); + + sizes.resize(nComp); + for (size_t i = 0; i < nComp; ++i) + sizes[i] = other.getNumRows(i); + + return *this; + } + /// Constructor for block-matrices BlockMapper(const SolverMatrix<Matrix<DOFMatrix* > >& sm ) : nComp(sm.getOriginalMat()->getSize()), @@ -106,6 +142,19 @@ namespace AMDiS { ncol = nrow; } + /// Constructor for block-matrices + BlockMapper(const Matrix<DOFMatrix*>& orMat ) + : nComp(orMat.getSize()), + rowOffset(0), colOffset(0), nrow(0), ncol(0), sizes(nComp) + { + const int ns = orMat.getNumRows(); + for (int i= 0; i < ns; i++) { + sizes[i] = orMat[i][i]->getFeSpace()->getAdmin()->getUsedSize(); + nrow += sizes[i]; + } + ncol = nrow; + } + /// Constructor for single matrix BlockMapper(const DOFMatrix* sm ) : nComp(1), diff --git a/AMDiS/src/solver/MatrixStreams.h b/AMDiS/src/solver/MatrixStreams.h index 0f7d62817573a0d2f1413643d7e5dc14b56d288e..9dc93235992224c5abe602aebaeb0dcc9879e5ad 100644 --- a/AMDiS/src/solver/MatrixStreams.h +++ b/AMDiS/src/solver/MatrixStreams.h @@ -34,7 +34,8 @@ namespace AMDiS { using namespace mtl::operations; - template< typename VecT, typename CurMap> + template< typename VecT, typename CurMap, + typename Updater = update_store<typename traits::category<VecT>::value_type> > struct VecMap { VecT& vec; CurMap& mapper; @@ -42,7 +43,7 @@ namespace AMDiS { vec(vec),mapper(mapper) {} }; - template< typename MatT, typename CurMap> + template< typename MatT, typename CurMap > struct MatMap { MatT& mat; CurMap& mapper; @@ -52,9 +53,9 @@ namespace AMDiS { // requires MatrixType to have an inserter template< typename MatrixType, typename Mapper > - void operator<<(MatrixType& matrix, MatMap<const SolverMatrix<Matrix<DOFMatrix* > >, Mapper >& Asolver) + void operator<<(MatrixType& matrix, MatMap<const Matrix<DOFMatrix*>, Mapper >& Asolver) { - const Matrix< DOFMatrix* >& A = *(Asolver.mat.getOriginalMat()); + const Matrix< DOFMatrix* >& A = Asolver.mat; int ns = A.getNumRows(); Mapper& mapper(Asolver.mapper); @@ -67,8 +68,8 @@ namespace AMDiS { nnz += A[rb][cb]->getBaseMatrix().nnz(); { - typedef mtl::matrix::mapped_inserter< typename Collection< MatrixType >::Inserter, Mapper > Inserter; - Inserter ins(matrix, mapper, int(1.2 * nnz / matrix.dim1())); + typedef typename Mapper::template Inserter< typename Collection<MatrixType>::Inserter >::type MappedInserter; + MappedInserter ins(matrix, mapper, int(1.2 * nnz / matrix.dim1())); for (int rb = 0; rb < ns; ++rb) { mapper.setRow(rb); for (int cb = 0; cb < ns; ++cb) { @@ -79,6 +80,14 @@ namespace AMDiS { } } } + + // requires MatrixType to have an inserter + template< typename MatrixType, typename Mapper > + void operator<<(MatrixType& matrix, MatMap<const SolverMatrix<Matrix<DOFMatrix* > >, Mapper >& Asolver) + { + MatMap<const Matrix<DOFMatrix* >, Mapper> mapMat(*Asolver.mat.getOriginalMat(), Asolver.mapper); + matrix << mapMat; + } // requires MatrixType to have an inserter @@ -98,8 +107,8 @@ namespace AMDiS { nnz += A[rb][cb]->nnz(); { - typedef mtl::matrix::mapped_inserter< typename Collection< MatrixType >::Inserter, Mapper > Inserter; - Inserter ins(matrix, mapper, int(1.2 * nnz / matrix.dim1())); + typedef typename Mapper::template Inserter< typename Collection<MatrixType>::Inserter >::type MappedInserter; + MappedInserter ins(matrix, mapper, int(1.2 * nnz / matrix.dim1())); for (int rb = 0; rb < ns; ++rb) { mapper.setRow(rb); for (int cb = 0; cb < ns; ++cb) { @@ -111,15 +120,16 @@ namespace AMDiS { } } - template< typename Vector, typename CurMap > - void operator<<(Vector& dest, VecMap<const DOFVector<double>, CurMap >& rhs) + template< typename Vector, typename CurMap, typename Updater > + void operator<<(Vector& dest, VecMap<const DOFVector<double>, CurMap, Updater>& rhs) { DOFConstIterator<double> it_x(&rhs.vec, USED_DOFS); - int counter(0); - typedef typename Collection< Vector >::Inserter Inserter ; + size_t counter(0); + typedef typename mtl::vector::inserter<Vector, Updater> Inserter; Inserter swapIns(dest); - mtl::vector::mapped_inserter< Inserter, CurMap > ins(swapIns, rhs.mapper); + typedef mtl::vector::mapped_inserter<Inserter, CurMap> MappedInserter; + MappedInserter ins(swapIns, rhs.mapper); for (it_x.reset(); !it_x.end(); ++it_x) { ins[counter] << *it_x; @@ -127,11 +137,11 @@ namespace AMDiS { } } - template< typename Vector, typename CurMap > - inline void operator>>(const Vector& dest, VecMap<DOFVector<double>, CurMap >& rhs) + template< typename Vector, typename CurMap, typename Updater > + inline void operator>>(const Vector& dest, VecMap<DOFVector<double>, CurMap, Updater>& rhs) { DOFVector<double>::Iterator it_x(&rhs.vec, USED_DOFS); - int counter(0); + size_t counter(0); { mtl::vector::extracter< Vector > extracter(dest); for (it_x.reset(); !it_x.end(); ++it_x) { @@ -141,18 +151,19 @@ namespace AMDiS { } } - template< typename Vector, typename CurMap > - void operator<<(Vector& dest, VecMap<const SystemVector, CurMap>& rhs) + template< typename Vector, typename CurMap, typename Updater > + void operator<<(Vector& dest, VecMap<const SystemVector, CurMap, Updater>& rhs) { int ns = rhs.vec.getSize(); // Number of systems. // Copy vectors - typedef typename Collection< Vector >::Inserter Inserter; + typedef typename mtl::vector::inserter<Vector, Updater> Inserter; Inserter swapInserter(dest); - mtl::vector::mapped_inserter< Inserter, CurMap > ins(swapInserter, rhs.mapper); + typedef mtl::vector::mapped_inserter<Inserter, CurMap> MappedInserter; + MappedInserter ins(swapInserter, rhs.mapper); for (int i = 0; i < ns; i++) { DOFConstIterator<double> it_source(rhs.vec.getDOFVector(i), USED_DOFS); - int counter(0); + size_t counter(0); rhs.mapper.setRow(i); for (it_source.reset(); !it_source.end(); ++it_source) { ins[counter] << *it_source; @@ -161,8 +172,8 @@ namespace AMDiS { } } - template< typename Vector, typename CurMap > - void operator>>(const Vector& dest, VecMap<SystemVector, CurMap>& rhs) + template< typename Vector, typename CurMap, typename Updater > + void operator>>(const Vector& dest, VecMap<SystemVector, CurMap, Updater>& rhs) { int ns = rhs.vec.getSize(); // Number of systems. diff --git a/AMDiS/src/utility/is_symmetric.hpp b/AMDiS/src/utility/is_symmetric.hpp index 8f030ccf85be33b9206c2f5e30ee6cd52a7eccc8..cccde03e351a1ac30fd03b4399f6a2318f67c760 100644 --- a/AMDiS/src/utility/is_symmetric.hpp +++ b/AMDiS/src/utility/is_symmetric.hpp @@ -33,12 +33,12 @@ namespace AMDiS { using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end; namespace traits= mtl::traits; - traits::row<Matrix>::type row(matrix); - traits::col<Matrix>::type col(matrix); - traits::const_value<Matrix>::type value(matrix); + typename traits::row<Matrix>::type row(matrix); + typename traits::col<Matrix>::type col(matrix); + typename traits::const_value<Matrix>::type value(matrix); - typedef traits::range_generator<major, Matrix>::type cursor_type; - typedef traits::range_generator<nz, cursor_type>::type icursor_type; + typedef typename traits::range_generator<major, Matrix>::type cursor_type; + typedef typename traits::range_generator<nz, cursor_type>::type icursor_type; for (cursor_type cursor = begin<major>(matrix), cend = end<major>(matrix); cursor != cend; ++cursor) for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor)