Commit 6577015f authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

* Bugfix and small changes to make it a little bit faster

parent 86fe284b
...@@ -24,7 +24,11 @@ namespace AMDiS { ...@@ -24,7 +24,11 @@ namespace AMDiS {
bool optimized, bool optimized,
FirstOrderType type) FirstOrderType type)
: SubAssembler(op, assembler, quad, 1, optimized, type) : SubAssembler(op, assembler, quad, 1, optimized, type)
{} {
VectorOfFixVecs<DimVec<double> >
Lb(assembler->getRowFESpace()->getMesh()->getDim(), 1, NO_INIT);
tmpLb.resize(omp_get_overall_max_threads(), Lb);
}
FirstOrderAssembler* FirstOrderAssembler*
FirstOrderAssembler::getSubAssembler(Operator* op, FirstOrderAssembler::getSubAssembler(Operator* op,
...@@ -95,8 +99,8 @@ namespace AMDiS { ...@@ -95,8 +99,8 @@ namespace AMDiS {
} else { } else {
newAssembler = newAssembler =
(type == GRD_PSI) ? (type == GRD_PSI) ?
dynamic_cast<FirstOrderAssembler*>( NEW Quad10(op, assembler, quad)) : dynamic_cast<FirstOrderAssembler*>(NEW Quad10(op, assembler, quad)) :
dynamic_cast<FirstOrderAssembler*>( NEW Quad01(op, assembler, quad)); dynamic_cast<FirstOrderAssembler*>(NEW Quad01(op, assembler, quad));
} }
} }
...@@ -108,23 +112,22 @@ namespace AMDiS { ...@@ -108,23 +112,22 @@ namespace AMDiS {
: FirstOrderAssembler(op, assembler, quad, false, GRD_PSI) : FirstOrderAssembler(op, assembler, quad, false, GRD_PSI)
{} {}
void Stand10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat) void Stand10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{ {
DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0); DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
double *phival = GET_MEMORY(double, nCol);
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts(); const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
const BasisFunction *phi = owner->getColFESpace()->getBasisFcts(); const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
int nPoints = quadrature->getNumPoints(); int nPoints = quadrature->getNumPoints();
int myRank = omp_get_thread_num();
VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
Lb.resize(nPoints);
std::vector<double> phival(nCol);
VectorOfFixVecs<DimVec<double> > Lb(dim, nPoints, NO_INIT);
for (int iq = 0; iq < nPoints; iq++) { for (int iq = 0; iq < nPoints; iq++) {
Lb[iq].set(0.0); Lb[iq].set(0.0);
} }
int myRank = omp_get_thread_num();
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) { for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb); (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
} }
...@@ -143,7 +146,6 @@ namespace AMDiS { ...@@ -143,7 +146,6 @@ namespace AMDiS {
} }
} }
} }
FREE_MEMORY(phival, double, nCol);
} }
void Stand10::calculateElementMatrix(const ElInfo *rowElInfo, void Stand10::calculateElementMatrix(const ElInfo *rowElInfo,
...@@ -165,13 +167,14 @@ namespace AMDiS { ...@@ -165,13 +167,14 @@ namespace AMDiS {
TEST_EXIT(m)("No subElemCoordsMat!\n"); TEST_EXIT(m)("No subElemCoordsMat!\n");
int nPoints = quadrature->getNumPoints(); int nPoints = quadrature->getNumPoints();
int myRank = omp_get_thread_num();
VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
Lb.resize(nPoints);
VectorOfFixVecs<DimVec<double> > Lb(dim, nPoints, NO_INIT);
for (int iq = 0; iq < nPoints; iq++) { for (int iq = 0; iq < nPoints; iq++) {
Lb[iq].set(0.0); Lb[iq].set(0.0);
} }
int myRank = omp_get_thread_num();
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) { for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<FirstOrderTerm*>((terms[myRank][i])))-> (static_cast<FirstOrderTerm*>((terms[myRank][i])))->
getLb(smallElInfo, nPoints, Lb); getLb(smallElInfo, nPoints, Lb);
...@@ -202,17 +205,41 @@ namespace AMDiS { ...@@ -202,17 +205,41 @@ namespace AMDiS {
} }
} }
void Stand10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
{
DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
int nPoints = quadrature->getNumPoints();
int myRank = omp_get_thread_num();
VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
Lb.resize(nPoints);
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq].set(0.0);
}
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
for (int i = 0; i < nRow; i++) {
(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
(*vec)[i] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi);
}
}
}
Quad10::Quad10(Operator *op, Assembler *assembler, Quadrature *quad) Quad10::Quad10(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, true, GRD_PSI) : FirstOrderAssembler(op, assembler, quad, true, GRD_PSI)
{ {}
}
void Quad10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat) void Quad10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{ {
VectorOfFixVecs<DimVec<double> > *grdPsi; VectorOfFixVecs<DimVec<double> > *grdPsi;
const double *phi;
if (firstCall) { if (firstCall) {
#ifdef _OPENMP #ifdef _OPENMP
...@@ -228,13 +255,14 @@ namespace AMDiS { ...@@ -228,13 +255,14 @@ namespace AMDiS {
} }
int nPoints = quadrature->getNumPoints(); int nPoints = quadrature->getNumPoints();
int myRank = omp_get_thread_num();
VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
Lb.resize(nPoints);
VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
for (int iq = 0; iq < nPoints; iq++) { for (int iq = 0; iq < nPoints; iq++) {
Lb[iq].set(0.0); Lb[iq].set(0.0);
} }
int myRank = omp_get_thread_num();
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) { for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb); (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
} }
...@@ -242,7 +270,7 @@ namespace AMDiS { ...@@ -242,7 +270,7 @@ namespace AMDiS {
for (int iq = 0; iq < nPoints; iq++) { for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet(); Lb[iq] *= elInfo->getDet();
phi = phiFast->getPhi(iq); const double *phi = phiFast->getPhi(iq);
grdPsi = psiFast->getGradient(iq); grdPsi = psiFast->getGradient(iq);
for (int i = 0; i < nRow; i++) { for (int i = 0; i < nRow; i++) {
...@@ -252,6 +280,45 @@ namespace AMDiS { ...@@ -252,6 +280,45 @@ namespace AMDiS {
} }
} }
void Quad10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
{
VectorOfFixVecs<DimVec<double> > *grdPsi;
if (firstCall) {
#ifdef _OPENMP
#pragma omp critical
#endif
{
const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
basFcts = owner->getColFESpace()->getBasisFcts();
phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
firstCall = false;
}
}
int nPoints = quadrature->getNumPoints();
int myRank = omp_get_thread_num();
VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
Lb.resize(nPoints);
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq].set(0.0);
}
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
grdPsi = psiFast->getGradient(iq);
for (int i = 0; i < nRow; i++) {
(*vec)[i] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]);
}
}
}
Pre10::Pre10(Operator *op, Assembler *assembler, Quadrature *quad) Pre10::Pre10(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, true, GRD_PSI) : FirstOrderAssembler(op, assembler, quad, true, GRD_PSI)
...@@ -261,7 +328,6 @@ namespace AMDiS { ...@@ -261,7 +328,6 @@ namespace AMDiS {
void Pre10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat) void Pre10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{ {
VectorOfFixVecs<DimVec<double> > Lb(dim, 1, NO_INIT);
const int *k; const int *k;
const double *values; const double *values;
...@@ -281,6 +347,8 @@ namespace AMDiS { ...@@ -281,6 +347,8 @@ namespace AMDiS {
const int **nEntries = q10->getNumberEntries(); const int **nEntries = q10->getNumberEntries();
int myRank = omp_get_thread_num(); int myRank = omp_get_thread_num();
// Do not need do resize Lb, because it's size is always at least one.
VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
Lb[0].set(0.0); Lb[0].set(0.0);
for (int i = 0; i < static_cast<int>( terms[myRank].size()); i++) { for (int i = 0; i < static_cast<int>( terms[myRank].size()); i++) {
...@@ -307,7 +375,6 @@ namespace AMDiS { ...@@ -307,7 +375,6 @@ namespace AMDiS {
: FirstOrderAssembler(op, assembler, quad, false, GRD_PHI) : FirstOrderAssembler(op, assembler, quad, false, GRD_PHI)
{} {}
void Stand01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat) void Stand01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{ {
VectorOfFixVecs<DimVec<double> > grdPhi(dim, nCol, NO_INIT); VectorOfFixVecs<DimVec<double> > grdPhi(dim, nCol, NO_INIT);
...@@ -316,8 +383,9 @@ namespace AMDiS { ...@@ -316,8 +383,9 @@ namespace AMDiS {
const BasisFunction *phi = owner->getColFESpace()->getBasisFcts(); const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
int nPoints = quadrature->getNumPoints(); int nPoints = quadrature->getNumPoints();
VectorOfFixVecs<DimVec<double> > Lb(dim, nPoints, NO_INIT);
int myRank = omp_get_thread_num(); int myRank = omp_get_thread_num();
VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
Lb.resize(nPoints);
for (int iq = 0; iq < nPoints; iq++) { for (int iq = 0; iq < nPoints; iq++) {
Lb[iq].set(0.0); Lb[iq].set(0.0);
...@@ -347,7 +415,7 @@ namespace AMDiS { ...@@ -347,7 +415,7 @@ namespace AMDiS {
const ElInfo *largeElInfo, const ElInfo *largeElInfo,
ElementMatrix *mat) ElementMatrix *mat)
{ {
FUNCNAME("Stand10::calculateElementMatrix()"); FUNCNAME("Stand01::calculateElementMatrix()");
TEST_EXIT((nRow <= 3) && (nCol <= 3))("not yet!\n"); TEST_EXIT((nRow <= 3) && (nCol <= 3))("not yet!\n");
...@@ -360,13 +428,14 @@ namespace AMDiS { ...@@ -360,13 +428,14 @@ namespace AMDiS {
TEST_EXIT(m)("No subElemCoordsMat!\n"); TEST_EXIT(m)("No subElemCoordsMat!\n");
int nPoints = quadrature->getNumPoints(); int nPoints = quadrature->getNumPoints();
int myRank = omp_get_thread_num();
VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
Lb.resize(nPoints);
VectorOfFixVecs<DimVec<double> > Lb(dim, nPoints, NO_INIT);
for (int iq = 0; iq < nPoints; iq++) { for (int iq = 0; iq < nPoints; iq++) {
Lb[iq].set(0.0); Lb[iq].set(0.0);
} }
int myRank = omp_get_thread_num();
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) { for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<FirstOrderTerm*>((terms[myRank][i])))-> (static_cast<FirstOrderTerm*>((terms[myRank][i])))->
getLb(smallElInfo, nPoints, Lb); getLb(smallElInfo, nPoints, Lb);
...@@ -396,35 +465,9 @@ namespace AMDiS { ...@@ -396,35 +465,9 @@ namespace AMDiS {
} }
} }
void Stand10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
{
DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
int nPoints = quadrature->getNumPoints();
VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
int myRank = omp_get_thread_num();
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq].set(0.0);
}
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
for (int i = 0; i < nRow; i++) {
(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
(*vec)[i] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi);
}
}
}
Quad01::Quad01(Operator *op, Assembler *assembler, Quadrature *quad) Quad01::Quad01(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, true, GRD_PHI) : FirstOrderAssembler(op, assembler, quad, true, GRD_PHI)
{ {}
}
void Quad01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat) void Quad01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{ {
...@@ -444,8 +487,9 @@ namespace AMDiS { ...@@ -444,8 +487,9 @@ namespace AMDiS {
} }
int nPoints = quadrature->getNumPoints(); int nPoints = quadrature->getNumPoints();
VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
int myRank = omp_get_thread_num(); int myRank = omp_get_thread_num();
VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
Lb.resize(nPoints);
for (int iq = 0; iq < nPoints; iq++) { for (int iq = 0; iq < nPoints; iq++) {
Lb[iq].set(0.0); Lb[iq].set(0.0);
...@@ -467,45 +511,6 @@ namespace AMDiS { ...@@ -467,45 +511,6 @@ namespace AMDiS {
} }
} }
void Quad10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
{
VectorOfFixVecs<DimVec<double> > *grdPsi;
if (firstCall) {
#ifdef _OPENMP
#pragma omp critical
#endif
{
const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
basFcts = owner->getColFESpace()->getBasisFcts();
phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
firstCall = false;
}
}
int nPoints = quadrature->getNumPoints();
VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
int myRank = omp_get_thread_num();
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq].set(0.0);
}
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
grdPsi = psiFast->getGradient(iq);
for (int i = 0; i < nRow; i++) {
(*vec)[i] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]);
}
}
}
Pre01::Pre01(Operator *op, Assembler *assembler, Quadrature *quad) Pre01::Pre01(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, true, GRD_PHI) : FirstOrderAssembler(op, assembler, quad, true, GRD_PHI)
{ {
...@@ -513,8 +518,6 @@ namespace AMDiS { ...@@ -513,8 +518,6 @@ namespace AMDiS {
void Pre01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat) void Pre01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{ {
VectorOfFixVecs<DimVec<double> > Lb(dim,1,NO_INIT);
const int *l; const int *l;
const double *values; const double *values;
...@@ -534,6 +537,8 @@ namespace AMDiS { ...@@ -534,6 +537,8 @@ namespace AMDiS {
const int **nEntries = q01->getNumberEntries(); const int **nEntries = q01->getNumberEntries();
int myRank = omp_get_thread_num(); int myRank = omp_get_thread_num();
// Do not need to resize Lb, because it's size is always at least one!
VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
Lb[0].set(0.0); Lb[0].set(0.0);
for (int i = 0; i < static_cast<int>( terms[myRank].size()); i++) { for (int i = 0; i < static_cast<int>( terms[myRank].size()); i++) {
...@@ -557,8 +562,6 @@ namespace AMDiS { ...@@ -557,8 +562,6 @@ namespace AMDiS {
void Pre10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec) void Pre10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
{ {
VectorOfFixVecs<DimVec<double> > Lb(dim,1,NO_INIT);
const int *k; const int *k;
const double *values; const double *values;
...@@ -578,6 +581,8 @@ namespace AMDiS { ...@@ -578,6 +581,8 @@ namespace AMDiS {
const int *nEntries = q1->getNumberEntries(); const int *nEntries = q1->getNumberEntries();
int myRank = omp_get_thread_num(); int myRank = omp_get_thread_num();
// Do not need to resize Lb, because it's size is always at least one!
VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
Lb[0].set(0.0); Lb[0].set(0.0);
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) { for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
......
...@@ -46,6 +46,12 @@ namespace AMDiS { ...@@ -46,6 +46,12 @@ namespace AMDiS {
protected: protected:
/** \brief
* Thread safe temporary vector of DimMats for calculation in
* function calculateElementMatrix().
*/
std::vector<VectorOfFixVecs<DimVec<double> > > tmpLb;
/// List of all yet created optimized zero order assemblers for grdPsi. /// List of all yet created optimized zero order assemblers for grdPsi.
static std::vector<SubAssembler*> optimizedSubAssemblersGrdPsi; static std::vector<SubAssembler*> optimizedSubAssemblersGrdPsi;
...@@ -71,7 +77,7 @@ namespace AMDiS { ...@@ -71,7 +77,7 @@ namespace AMDiS {
public: public:
MEMORY_MANAGED(Stand10); MEMORY_MANAGED(Stand10);
/// Constructor. /// Constructor
Stand10(Operator *op, Assembler *assembler, Quadrature *quad); Stand10(Operator *op, Assembler *assembler, Quadrature *quad);
/// Implements SubAssembler::calculateElementMatrix(). /// Implements SubAssembler::calculateElementMatrix().
......
...@@ -108,33 +108,23 @@ namespace AMDiS { ...@@ -108,33 +108,23 @@ namespace AMDiS {
this->set(ini); this->set(ini);
} }
/** \brief /// Initialisation for dim.
* Initialisation for dim. inline void init(int dim) {
*/
inline void init(int dim)
{
this->resize(calcSize(dim)); this->resize(calcSize(dim));
} }
/** \brief /// Initialisation for size
* Initialisation for size inline void initSize(int size) {
*/
inline void initSize(int size)
{
this->resize(size); this->resize(size);
} }
/** \brief /// Returns the \ref size_ of the FixVec.
* Returns the \ref size_ of the FixVec.
*/
inline int size() const { inline int size() const {
return this->getSize(); return this->getSize();
} }
protected: protected:
/** \brief /// Determines needed vector size.
* Determines needed vector size.
*/
static int calcSize(int dim) { static int calcSize(int dim) {
if (dim < 0) { if (dim < 0) {
return Global::getGeo(WORLD); return Global::getGeo(WORLD);
...@@ -172,15 +162,15 @@ namespace AMDiS { ...@@ -172,15 +162,15 @@ namespace AMDiS {
* FixVec's constructors. size_ is the number of contained FixVecs. initType * FixVec's constructors. size_ is the number of contained FixVecs. initType