Commit ba2f41ec authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

And more performance optimizations.

parent 2d7b92be
...@@ -243,7 +243,6 @@ namespace AMDiS { ...@@ -243,7 +243,6 @@ namespace AMDiS {
("invalid basis functions"); ("invalid basis functions");
int dow = Global::getGeo(WORLD); int dow = Global::getGeo(WORLD);
int myRank = omp_get_thread_num();
int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints(); int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
static WorldVector<double> *grd = NULL; static WorldVector<double> *grd = NULL;
WorldVector<double> *result; WorldVector<double> *result;
...@@ -261,8 +260,8 @@ namespace AMDiS { ...@@ -261,8 +260,8 @@ namespace AMDiS {
result = grd; result = grd;
} }
const ElementVector& localVec = localVectors[myRank]; ElementVector localVec(nBasFcts);
getLocalVector(elInfo->getElement(), const_cast<ElementVector&>(localVec)); getLocalVector(elInfo->getElement(), localVec);
mtl::dense_vector<double> grd1(dim + 1); mtl::dense_vector<double> grd1(dim + 1);
int parts = Global::getGeo(PARTS, dim); int parts = Global::getGeo(PARTS, dim);
...@@ -327,7 +326,6 @@ namespace AMDiS { ...@@ -327,7 +326,6 @@ namespace AMDiS {
("invalid basis functions"); ("invalid basis functions");
int dow = Global::getGeo(WORLD); int dow = Global::getGeo(WORLD);
int myRank = omp_get_thread_num();
int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints(); int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
static WorldVector<double> *grd = NULL; static WorldVector<double> *grd = NULL;
WorldVector<double> *result; WorldVector<double> *result;
...@@ -345,8 +343,8 @@ namespace AMDiS { ...@@ -345,8 +343,8 @@ namespace AMDiS {
result = grd; result = grd;
} }
const ElementVector& localVec = localVectors[myRank]; ElementVector localVec(nBasFcts);
getLocalVector(largeElInfo->getElement(), const_cast<ElementVector&>(localVec)); getLocalVector(largeElInfo->getElement(), localVec);
const BasisFunction *basFcts = feSpace->getBasisFcts(); const BasisFunction *basFcts = feSpace->getBasisFcts();
mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree()); mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
...@@ -405,7 +403,6 @@ namespace AMDiS { ...@@ -405,7 +403,6 @@ namespace AMDiS {
Element *el = elInfo->getElement(); Element *el = elInfo->getElement();
int myRank = omp_get_thread_num();
int dow = Global::getGeo(WORLD); int dow = Global::getGeo(WORLD);
int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints(); int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
...@@ -423,8 +420,8 @@ namespace AMDiS { ...@@ -423,8 +420,8 @@ namespace AMDiS {
result = vec; result = vec;
} }
const ElementVector& localVec = localVectors[myRank]; ElementVector localVec(nBasFcts);
getLocalVector(el, const_cast<ElementVector&>(localVec)); getLocalVector(el, localVec);
DimMat<double> D2Tmp(dim, DEFAULT_VALUE, 0.0); DimMat<double> D2Tmp(dim, DEFAULT_VALUE, 0.0);
int parts = Global::getGeo(PARTS, dim); int parts = Global::getGeo(PARTS, dim);
...@@ -452,7 +449,7 @@ namespace AMDiS { ...@@ -452,7 +449,7 @@ namespace AMDiS {
} }
} else { } else {
const BasisFunction *basFcts = feSpace->getBasisFcts(); const BasisFunction *basFcts = feSpace->getBasisFcts();
DimMat<double>* D2Phi = D2Phis[omp_get_thread_num()]; DimMat<double> D2Phi(dim, NO_INIT);
for (int iq = 0; iq < nPoints; iq++) { for (int iq = 0; iq < nPoints; iq++) {
for (int k = 0; k < parts; k++) for (int k = 0; k < parts; k++)
...@@ -461,11 +458,11 @@ namespace AMDiS { ...@@ -461,11 +458,11 @@ namespace AMDiS {
for (int i = 0; i < nBasFcts; i++) { for (int i = 0; i < nBasFcts; i++) {
WARNING("not tested after index correction\n"); WARNING("not tested after index correction\n");
(*(basFcts->getD2Phi(i)))(quad->getLambda(iq), *D2Phi); (*(basFcts->getD2Phi(i)))(quad->getLambda(iq), D2Phi);
for (int k = 0; k < parts; k++) for (int k = 0; k < parts; k++)
for (int l = 0; l < parts; l++) for (int l = 0; l < parts; l++)
D2Tmp[k][l] += localVec[i] * (*D2Phi)[k][l]; D2Tmp[k][l] += localVec[i] * D2Phi[k][l];
} }
for (int i = 0; i < dow; i++) for (int i = 0; i < dow; i++)
...@@ -496,7 +493,7 @@ namespace AMDiS { ...@@ -496,7 +493,7 @@ namespace AMDiS {
this->set(0.0); this->set(0.0);
DegreeOfFreedom *myLocalIndices = localIndices[omp_get_thread_num()]; std::vector<DegreeOfFreedom> myLocalIndices(nBasisFcts);
ElementVector sourceLocalCoeffs(nSourceBasisFcts); ElementVector sourceLocalCoeffs(nSourceBasisFcts);
if (feSpace->getMesh() == sourceFeSpace->getMesh()) { if (feSpace->getMesh() == sourceFeSpace->getMesh()) {
...@@ -582,11 +579,11 @@ namespace AMDiS { ...@@ -582,11 +579,11 @@ namespace AMDiS {
const BasisFunction *basFcts = feSpace->getBasisFcts(); const BasisFunction *basFcts = feSpace->getBasisFcts();
const BasisFunction *vBasFcts = vFeSpace->getBasisFcts(); const BasisFunction *vBasFcts = vFeSpace->getBasisFcts();
int numBasFcts = basFcts->getNumber(); int nBasFcts = basFcts->getNumber();
int vNumBasFcts = vBasFcts->getNumber(); int vNumBasFcts = vBasFcts->getNumber();
if (feSpace->getMesh() == vFeSpace->getMesh()) { if (feSpace->getMesh() == vFeSpace->getMesh()) {
DegreeOfFreedom *myLocalIndices = localIndices[omp_get_thread_num()]; std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
mtl::dense_vector<WorldVector<double> > vLocalCoeffs(vNumBasFcts); mtl::dense_vector<WorldVector<double> > vLocalCoeffs(vNumBasFcts);
Mesh *mesh = feSpace->getMesh(); Mesh *mesh = feSpace->getMesh();
TraverseStack stack; TraverseStack stack;
...@@ -598,7 +595,7 @@ namespace AMDiS { ...@@ -598,7 +595,7 @@ namespace AMDiS {
basFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices); basFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
v->getLocalVector(el, vLocalCoeffs); v->getLocalVector(el, vLocalCoeffs);
for (int i = 0; i < numBasFcts; i++) { for (int i = 0; i < nBasFcts; i++) {
if (vec[myLocalIndices[i]] == nul) { if (vec[myLocalIndices[i]] == nul) {
coords = basFcts->getCoords(i); coords = basFcts->getCoords(i);
vec[myLocalIndices[i]] += vBasFcts->evalUh(*coords, vLocalCoeffs, NULL) * factor; vec[myLocalIndices[i]] += vBasFcts->evalUh(*coords, vLocalCoeffs, NULL) * factor;
...@@ -776,9 +773,8 @@ namespace AMDiS { ...@@ -776,9 +773,8 @@ namespace AMDiS {
vecAtQPs.change_dim(nPoints); vecAtQPs.change_dim(nPoints);
const ElementVector &localVec = localVectors[omp_get_thread_num()]; ElementVector localVec(nBasFcts);
getLocalVector(largeElInfo->getElement(), const_cast<ElementVector&>(localVec)); getLocalVector(largeElInfo->getElement(), localVec);
mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree()); mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
for (int i = 0; i < nPoints; i++) { for (int i = 0; i < nPoints; i++) {
......
...@@ -64,9 +64,6 @@ namespace AMDiS { ...@@ -64,9 +64,6 @@ namespace AMDiS {
virtual ~DOFVectorBase(); virtual ~DOFVectorBase();
/// Creates temporary used data structures.
void createTempData();
/** \brief /** \brief
* For the given element, this function returns an array of all DOFs of this * For the given element, this function returns an array of all DOFs of this
* DOFVector that are defined on this element. * DOFVector that are defined on this element.
...@@ -263,15 +260,6 @@ namespace AMDiS { ...@@ -263,15 +260,6 @@ namespace AMDiS {
/// Number of basis functions of the used finite element space. /// Number of basis functions of the used finite element space.
int nBasFcts; int nBasFcts;
/// Are used to store temporary local dofs of an element. Thread safe.
std::vector<DegreeOfFreedom*> localIndices;
/// Are used to store temporary local values of an element. Thread safe.
std::vector<mtl::dense_vector<T> > localVectors;
/// Temporarly used in \ref getD2AtQPs. Thread safe.
std::vector<DimMat<double>*> D2Phis;
/// Dimension of the mesh this DOFVectorBase belongs to /// Dimension of the mesh this DOFVectorBase belongs to
int dim; int dim;
...@@ -359,8 +347,6 @@ namespace AMDiS { ...@@ -359,8 +347,6 @@ namespace AMDiS {
this->name = rhs.name + "copy"; this->name = rhs.name + "copy";
if (this->feSpace && this->feSpace->getAdmin()) if (this->feSpace && this->feSpace->getAdmin())
(dynamic_cast<DOFAdmin*>(this->feSpace->getAdmin()))->addDOFIndexed(this); (dynamic_cast<DOFAdmin*>(this->feSpace->getAdmin()))->addDOFIndexed(this);
this->createTempData();
} }
/// Destructor /// Destructor
......
...@@ -31,7 +31,6 @@ ...@@ -31,7 +31,6 @@
#include "AbstractFunction.h" #include "AbstractFunction.h"
#include "BoundaryManager.h" #include "BoundaryManager.h"
#include "Assembler.h" #include "Assembler.h"
#include "OpenMP.h"
#include "Operator.h" #include "Operator.h"
#include "Parameters.h" #include "Parameters.h"
#include "Traverse.h" #include "Traverse.h"
...@@ -88,36 +87,14 @@ namespace AMDiS { ...@@ -88,36 +87,14 @@ namespace AMDiS {
{ {
nBasFcts = feSpace->getBasisFcts()->getNumber(); nBasFcts = feSpace->getBasisFcts()->getNumber();
dim = feSpace->getMesh()->getDim(); dim = feSpace->getMesh()->getDim();
createTempData();
} }
template<typename T> template<typename T>
DOFVectorBase<T>::~DOFVectorBase() DOFVectorBase<T>::~DOFVectorBase()
{ {}
for (int i = 0; i < static_cast<int>(localIndices.size()); i++) {
delete [] localIndices[i];
delete D2Phis[i];
}
}
template<typename T>
void DOFVectorBase<T>::createTempData()
{
localIndices.resize(omp_get_overall_max_threads());
localVectors.resize(omp_get_overall_max_threads());
D2Phis.resize(omp_get_overall_max_threads());
for (int i = 0; i < omp_get_overall_max_threads(); i++) {
localIndices[i] = new DegreeOfFreedom[this->nBasFcts];
localVectors[i].change_dim(this->nBasFcts);
D2Phis[i] = new DimMat<double>(dim, NO_INIT);
}
}
template<typename T> template<typename T>
DOFVector<T>::DOFVector(const FiniteElemSpace* f, std::string n) DOFVector<T>::DOFVector(const FiniteElemSpace* f, std::string n)
: DOFVectorBase<T>(f, n), : DOFVectorBase<T>(f, n),
...@@ -465,9 +442,10 @@ namespace AMDiS { ...@@ -465,9 +442,10 @@ namespace AMDiS {
const_cast<BasisFunction*>(basFct)->interpol(elinfo, 0, NULL, const_cast<BasisFunction*>(basFct)->interpol(elinfo, 0, NULL,
traverseVector->interFct, NULL); traverseVector->interFct, NULL);
DegreeOfFreedom *myLocalIndices = this->localIndices[omp_get_thread_num()];
basFct->getLocalIndices(const_cast<Element*>(elinfo->getElement()), admin, myLocalIndices);
int nBasFcts = basFct->getNumber(); int nBasFcts = basFct->getNumber();
std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
basFct->getLocalIndices(const_cast<Element*>(elinfo->getElement()),
admin, myLocalIndices);
for (int i = 0; i < nBasFcts; i++) for (int i = 0; i < nBasFcts; i++)
(*traverseVector)[myLocalIndices[i]] = inter_val[i]; (*traverseVector)[myLocalIndices[i]] = inter_val[i];
} }
...@@ -991,19 +969,14 @@ namespace AMDiS { ...@@ -991,19 +969,14 @@ namespace AMDiS {
("Element is defined on a different mesh than the DOF vector!\n"); ("Element is defined on a different mesh than the DOF vector!\n");
const DOFAdmin* admin = feSpace->getAdmin(); const DOFAdmin* admin = feSpace->getAdmin();
#ifdef _OPENMP
std::vector<DegreeOfFreedom> localIndices(nBasFcts);
feSpace->getBasisFcts()->getLocalIndices(el, admin, &(localIndices[0]));
#else
const DegreeOfFreedom *localIndices = const DegreeOfFreedom *localIndices =
feSpace->getBasisFcts()->getLocalIndices(el, admin, NULL); feSpace->getBasisFcts()->getLocalIndices(el, admin, NULL);
#endif
for (int i = 0; i < nBasFcts; i++) for (int i = 0; i < nBasFcts; i++)
d[i] = (*this)[localIndices[i]]; d[i] = (*this)[localIndices[i]];
} }
template<typename T> template<typename T>
void DOFVectorBase<T>::getVecAtQPs(const ElInfo *elInfo, void DOFVectorBase<T>::getVecAtQPs(const ElInfo *elInfo,
const Quadrature *quad, const Quadrature *quad,
...@@ -1019,9 +992,10 @@ namespace AMDiS { ...@@ -1019,9 +992,10 @@ namespace AMDiS {
TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts()) TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
("Invalid basis functions!"); ("Invalid basis functions!");
const mtl::dense_vector<T>& localVec = localVectors[omp_get_thread_num()]; const BasisFunction *basFcts = feSpace->getBasisFcts();
getLocalVector(elInfo->getElement(), int nBasFcts = basFcts->getNumber();
const_cast<mtl::dense_vector<T>& >(localVec)); mtl::dense_vector<T> localVec(nBasFcts);
getLocalVector(elInfo->getElement(), localVec);
if (quadFast) { if (quadFast) {
const mtl::dense2D<double>& phi = quadFast->getPhi(); const mtl::dense2D<double>& phi = quadFast->getPhi();
...@@ -1029,9 +1003,7 @@ namespace AMDiS { ...@@ -1029,9 +1003,7 @@ namespace AMDiS {
vecAtQPs = phi * localVec; vecAtQPs = phi * localVec;
} else { } else {
const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad; const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad;
const BasisFunction *basFcts = feSpace->getBasisFcts();
int nPoints = quadrature->getNumPoints(); int nPoints = quadrature->getNumPoints();
int nBasFcts = basFcts->getNumber();
vecAtQPs.change_dim(nPoints); vecAtQPs.change_dim(nPoints);
for (int i = 0; i < nPoints; i++) { for (int i = 0; i < nPoints; i++) {
......
...@@ -207,7 +207,7 @@ namespace AMDiS { ...@@ -207,7 +207,7 @@ namespace AMDiS {
for (unsigned int i = 0; i < terms.size(); i++) for (unsigned int i = 0; i < terms.size(); i++)
(static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, Lb); (static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, Lb);
const mtl::dense2D<double>& phi(phiFast->getPhi()); const mtl::dense2D<double>& phi = phiFast->getPhi();
for (int iq = 0; iq < nPoints; iq++) { for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet(); Lb[iq] *= elInfo->getDet();
...@@ -220,8 +220,6 @@ namespace AMDiS { ...@@ -220,8 +220,6 @@ namespace AMDiS {
mat[i][j] += factor * phi[iq][j] * dot(Lb[iq], grdPsi[i]); mat[i][j] += factor * phi[iq][j] * dot(Lb[iq], grdPsi[i]);
} }
} }
// TODO: Do the same performance update es in Quad01::calculateElementMatrix
} }
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment