Skip to content

GitLab

Commit 6eda4c4f authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

* Assembler now faster for dendrite code

parent f615ae6e
Hide whitespace changes
Inline Side-by-side
Showing with 335 additions and 344 deletions
AMDiS/src/Assembler.cc
View file @ 6eda4c4f
......@@ -122,7 +122,7 @@ namespace AMDiS {
{
Quadrature *localQuad = quad ? quad : quadrature;
const int numPoints = localQuad->getNumPoints();
const int nPoints = localQuad->getNumPoints();
// already calculated for this element ?
if (coordsValid) {
......@@ -130,19 +130,19 @@ namespace AMDiS {
}
if (coordsAtQPs) {
if (coordsNumAllocated != numPoints) {
if (coordsNumAllocated != nPoints) {
DELETE [] coordsAtQPs;
coordsAtQPs = NEW WorldVector<double>[numPoints];
coordsNumAllocated = numPoints;
coordsAtQPs = NEW WorldVector<double>[nPoints];
coordsNumAllocated = nPoints;
}
} else {
coordsAtQPs = NEW WorldVector<double>[numPoints];
coordsNumAllocated = numPoints;
coordsAtQPs = NEW WorldVector<double>[nPoints];
coordsNumAllocated = nPoints;
}
// set new values
WorldVector<double>* k = &(coordsAtQPs[0]);
for (int l = 0; k < &(coordsAtQPs[numPoints]); ++k, ++l) {
for (int l = 0; k < &(coordsAtQPs[nPoints]); ++k, ++l) {
elInfo->coordToWorld(localQuad->getLambda(l), k);
}
......@@ -223,7 +223,7 @@ namespace AMDiS {
Quadrature *localQuad = quad ? quad : quadrature;
if(!gradientsAtQPs[vec]) {
if (!gradientsAtQPs[vec]) {
gradientsAtQPs[vec] = new GradientsAtQPs;
}
gradientsAtQPs[vec]->values.resize(localQuad->getNumPoints());
......@@ -249,7 +249,7 @@ namespace AMDiS {
const BasisFunction *basFcts = vec->getFESpace()->getBasisFcts();
if (opt && !quad && sameFESpaces) {
if(psiFast->getBasisFunctions() == basFcts) {
if (psiFast->getBasisFunctions() == basFcts) {
vec->getGrdAtQPs(elInfo, NULL, psiFast, values);
} else {
vec->getGrdAtQPs(elInfo, NULL, phiFast, values);
......@@ -500,21 +500,21 @@ namespace AMDiS {
double psival;
double *phival = GET_MEMORY(double, nCol);
int numPoints = quadrature->getNumPoints();
double *c = GET_MEMORY(double, numPoints);
int nPoints = quadrature->getNumPoints();
double *c = GET_MEMORY(double, nPoints);
for (int iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++) {
c[iq] = 0.0;
}
int myRank = omp_get_thread_num();
::std::vector<OperatorTerm*>::iterator termIt;
for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, numPoints, c);
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
}
if (symmetric) {
for (int iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++) {
c[iq] *= elInfo->getDet();
// calculate phi at QPs only once!
......@@ -533,7 +533,7 @@ namespace AMDiS {
}
}
} else { // non symmetric assembling
for (int iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++) {
c[iq] *= elInfo->getDet();
// calculate phi at QPs only once!
......@@ -551,26 +551,26 @@ namespace AMDiS {
}
FREE_MEMORY(phival, double, nCol);
FREE_MEMORY(c, double, numPoints);
FREE_MEMORY(c, double, nPoints);
}
void Stand0::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
{
int numPoints = quadrature->getNumPoints();
int nPoints = quadrature->getNumPoints();
double *c = GET_MEMORY(double, numPoints);
double *c = GET_MEMORY(double, nPoints);
for (int iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++) {
c[iq] = 0.0;
}
int myRank = omp_get_thread_num();
::std::vector<OperatorTerm*>::iterator termIt;
for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, numPoints, c);
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
}
for (int iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++) {
c[iq] *= elInfo->getDet();
for (int i = 0; i < nRow; i++) {
......@@ -580,19 +580,29 @@ namespace AMDiS {
}
}
FREE_MEMORY(c, double, numPoints);
FREE_MEMORY(c, double, nPoints);
}
Quad0::Quad0(Operator *op, Assembler *assembler, Quadrature *quad)
: ZeroOrderAssembler(op, assembler, quad, true)
{
cPtrs.resize(omp_get_max_threads());
}
Quad0::~Quad0()
{
for (int i = 0; i < omp_get_max_threads(); i++) {
FREE_MEMORY(cPtrs[i], double, quadrature->getNumPoints());
}
}
void Quad0::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
const double *psi, *phi;
int nPoints = quadrature->getNumPoints();
int myRank = omp_get_thread_num();
if (firstCall) {
cPtrs[myRank] = GET_MEMORY(double, nPoints);
const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
basFcts = owner->getColFESpace()->getBasisFcts();
......@@ -600,25 +610,22 @@ namespace AMDiS {
firstCall = false;
}
int numPoints = quadrature->getNumPoints();
double *c = GET_MEMORY(double, numPoints);
for (int iq = 0; iq < numPoints; iq++) {
double *c = cPtrs[myRank];
for (int iq = 0; iq < nPoints; iq++) {
c[iq] = 0.0;
}
int myRank = omp_get_thread_num();
::std::vector<OperatorTerm*>::iterator termIt;
for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, numPoints, c);
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
}
if (symmetric) {
for (int iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++) {
c[iq] *= elInfo->getDet();
psi = psiFast->getPhi(iq);
phi = phiFast->getPhi(iq);
const double *psi = psiFast->getPhi(iq);
const double *phi = phiFast->getPhi(iq);
for (int i = 0; i < nRow; i++) {
(*mat)[i][i] += quadrature->getWeight(iq) * c[iq] * psi[i] * phi[i];
for (int j = i + 1; j < nCol; j++) {
......@@ -629,11 +636,11 @@ namespace AMDiS {
}
}
} else { /* non symmetric assembling */
for (int iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++) {
c[iq] *= elInfo->getDet();
psi = psiFast->getPhi(iq);
phi = phiFast->getPhi(iq);
const double *psi = psiFast->getPhi(iq);
const double *phi = phiFast->getPhi(iq);
for (int i = 0; i < nRow; i++) {
for (int j = 0; j < nCol; j++) {
(*mat)[i][j] += quadrature->getWeight(iq) * c[iq] * psi[i] * phi[j];
......@@ -641,7 +648,6 @@ namespace AMDiS {
}
}
}
FREE_MEMORY(c, double, numPoints);
}
void Quad0::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
......@@ -654,20 +660,20 @@ namespace AMDiS {
firstCall = false;
}
int numPoints = quadrature->getNumPoints();
double *c = GET_MEMORY(double, numPoints);
int nPoints = quadrature->getNumPoints();
double *c = GET_MEMORY(double, nPoints);
for (int iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++) {
c[iq] = 0.0;
}
int myRank = omp_get_thread_num();
::std::vector<OperatorTerm*>::iterator termIt;
for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, numPoints, c);
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
}
for (int iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++) {
c[iq] *= elInfo->getDet();
const double *psi = psiFast->getPhi(iq);
......@@ -675,7 +681,7 @@ namespace AMDiS {
(*vec)[i] += quadrature->getWeight(iq) * c[iq] * psi[i];
}
}
FREE_MEMORY(c, double, numPoints);
FREE_MEMORY(c, double, nPoints);
}
Pre0::Pre0(Operator *op, Assembler *assembler, Quadrature *quad)
......@@ -759,19 +765,19 @@ namespace AMDiS {
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
int numPoints = quadrature->getNumPoints();
int nPoints = quadrature->getNumPoints();
VectorOfFixVecs<DimVec<double> > Lb(dim, numPoints, NO_INIT);
for (int iq = 0; iq < numPoints; iq++) {
VectorOfFixVecs<DimVec<double> > Lb(dim, nPoints, NO_INIT);
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq].set(0.0);
}
int myRank = omp_get_thread_num();
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, numPoints, Lb);
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
for (int iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
for (int i = 0; i < nCol; i++) {
......@@ -808,19 +814,19 @@ namespace AMDiS {
firstCall = false;
}
int numPoints = quadrature->getNumPoints();
int nPoints = quadrature->getNumPoints();
VectorOfFixVecs<DimVec<double> > Lb(dim,numPoints,NO_INIT);
for (int iq = 0; iq < numPoints; iq++) {
VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq].set(0.0);
}
int myRank = omp_get_thread_num();
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, numPoints, Lb);
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
for (int iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
phi = phiFast->getPhi(iq);
......@@ -891,18 +897,18 @@ namespace AMDiS {
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
int numPoints = quadrature->getNumPoints();
VectorOfFixVecs<DimVec<double> > Lb(dim, numPoints, NO_INIT);
int nPoints = quadrature->getNumPoints();
VectorOfFixVecs<DimVec<double> > Lb(dim, nPoints, NO_INIT);
int myRank = omp_get_thread_num();
for (int iq = 0; iq < numPoints; iq++) {
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, numPoints, Lb);
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
for (int iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
for (int i = 0; i < nCol; i++) {
......@@ -921,18 +927,18 @@ namespace AMDiS {
{
DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
int numPoints = quadrature->getNumPoints();
VectorOfFixVecs<DimVec<double> > Lb(dim,numPoints,NO_INIT);
int nPoints = quadrature->getNumPoints();
VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
int myRank = omp_get_thread_num();
for (int iq = 0; iq < numPoints; iq++) {
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, numPoints, Lb);
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
for (int iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
for (int i = 0; i < nRow; i++) {
......@@ -959,18 +965,18 @@ namespace AMDiS {
firstCall = false;
}
int numPoints = quadrature->getNumPoints();
VectorOfFixVecs<DimVec<double> > Lb(dim,numPoints,NO_INIT);
int nPoints = quadrature->getNumPoints();
VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
int myRank = omp_get_thread_num();
for (int iq = 0; iq < numPoints; iq++) {
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, numPoints, Lb);
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
for (int iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
const double *psi = psiFast->getPhi(iq);
......@@ -995,18 +1001,18 @@ namespace AMDiS {
firstCall = false;
}
int numPoints = quadrature->getNumPoints();
VectorOfFixVecs<DimVec<double> > Lb(dim,numPoints,NO_INIT);
int nPoints = quadrature->getNumPoints();
VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
int myRank = omp_get_thread_num();
for (int iq = 0; iq < numPoints; iq++) {
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, numPoints, Lb);
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
for (int iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
grdPsi = psiFast->getGradient(iq);
......@@ -1179,14 +1185,32 @@ namespace AMDiS {
Quad2::Quad2(Operator *op, Assembler *assembler, Quadrature *quad)
: SecondOrderAssembler(op, assembler, quad, true)
{}
{
tmpLALt.resize(omp_get_max_threads());
}
Quad2::~Quad2()
{
int nPoints = quadrature->getNumPoints();
for (int i = 0; i < static_cast<int>(tmpLALt.size()); i++) {
for (int j = 0; j < nPoints; j++) {
DELETE tmpLALt[i][j];
}
DELETE [] tmpLALt[i];
}
}
void Quad2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
double val;
VectorOfFixVecs<DimVec<double> > *grdPsi, *grdPhi;
int nPoints = quadrature->getNumPoints();
int myRank = omp_get_thread_num();
if (firstCall) {
tmpLALt[myRank] = NEW DimMat<double>*[nPoints];
for (int j = 0; j < nPoints; j++) {
tmpLALt[myRank][j] = NEW DimMat<double>(dim, NO_INIT);
}
const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
basFcts = owner->getColFESpace()->getBasisFcts();
......@@ -1194,20 +1218,18 @@ namespace AMDiS {
firstCall = false;
}
int nPoints = quadrature->getNumPoints();
DimMat<double> **LALt = NEW DimMat<double>*[nPoints];
int myRank = omp_get_thread_num();
DimMat<double> **LALt = tmpLALt[myRank];
for (int i = 0; i < nPoints;i++) {
LALt[i] = NEW DimMat<double>(dim, NO_INIT);
}
for (int iq = 0; iq < nPoints; iq++) {
LALt[iq]->set(0.0);
for (int i = 0; i < nPoints; i++) {
LALt[i]->set(0.0);
}
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, nPoints, LALt);
}
VectorOfFixVecs<DimVec<double> > *grdPsi, *grdPhi;
if (symmetric) {
for (int iq = 0; iq < nPoints; iq++) {
(*LALt[iq]) *= elInfo->getDet();
......@@ -1220,14 +1242,13 @@ namespace AMDiS {
((*grdPsi)[i] * ((*LALt[iq]) * (*grdPhi)[i]));
for (int j = i + 1; j < nCol; j++) {
val = quadrature->getWeight(iq) * ((*grdPsi)[i] * ((*LALt[iq]) * (*grdPhi)[j]));
double val = quadrature->getWeight(iq) * ((*grdPsi)[i] * ((*LALt[iq]) * (*grdPhi)[j]));
(*mat)[i][j] += val;
(*mat)[j][i] += val;
}
}
}
}
else { /* non symmetric assembling */
} else { /* non symmetric assembling */
for (int iq = 0; iq < nPoints; iq++) {
(*LALt[iq]) *= elInfo->getDet();
......@@ -1242,11 +1263,6 @@ namespace AMDiS {
}
}
}
for (int i = 0;i < nPoints; i++)
DELETE LALt[i];
DELETE [] LALt;
}
Stand2::Stand2(Operator *op, Assembler *assembler, Quadrature *quad)
......
AMDiS/src/Assembler.h
View file @ 6eda4c4f
......@@ -117,7 +117,7 @@ namespace AMDiS {
/** \brief
* Returns \ref terms
*/
inline ::std::vector<OperatorTerm*> *getTerms() {
inline std::vector<OperatorTerm*> *getTerms() {
return &terms[omp_get_thread_num()];
};
......@@ -227,12 +227,12 @@ namespace AMDiS {
/** \brief
* Used for \ref getVectorAtQPs().
*/
::std::map<const DOFVectorBase<double>*, ValuesAtQPs*> valuesAtQPs;
std::map<const DOFVectorBase<double>*, ValuesAtQPs*> valuesAtQPs;
/** \brief
* Used for \ref getGradientsAtQPs().
*/
::std::map<const DOFVectorBase<double>*, GradientsAtQPs*> gradientsAtQPs;
std::map<const DOFVectorBase<double>*, GradientsAtQPs*> gradientsAtQPs;
/** \brief
* Set and updated by \ref initElement() for each ElInfo.
......@@ -278,7 +278,7 @@ namespace AMDiS {
/** \brief
* List of all terms with a contribution to this SubAssembler
*/
::std::vector< ::std::vector<OperatorTerm*> > terms;
std::vector< std::vector<OperatorTerm*> > terms;
/** \brief
*
......@@ -339,12 +339,12 @@ namespace AMDiS {
/** \brief
* List of all yet created optimized SubAssembler objects.
*/
static ::std::vector<SubAssembler*> optimizedSubAssemblers;
static std::vector<SubAssembler*> optimizedSubAssemblers;
/** \brief
* List of all yet created standard SubAssembler objects.
*/
static ::std::vector<SubAssembler*> standardSubAssemblers;
static std::vector<SubAssembler*> standardSubAssemblers;
};
// ============================================================================
......@@ -398,6 +398,8 @@ namespace AMDiS {
*/
Quad0(Operator *op, Assembler *assembler, Quadrature *quad);
~Quad0();
/** \brief
* Implements SubAssembler::calculateElementMatrix().
*/
......@@ -407,6 +409,9 @@ namespace AMDiS {
* Implements SubAssembler::calculateElementVector().
*/
void calculateElementVector(const ElInfo *elInfo, ElementVector *vec);
protected:
std::vector<double*> cPtrs;
};
......@@ -502,22 +507,22 @@ namespace AMDiS {
/** \brief
* List of all yet created optimized zero order assemblers for grdPsi.
*/
static ::std::vector<SubAssembler*> optimizedSubAssemblersGrdPsi;
static std::vector<SubAssembler*> optimizedSubAssemblersGrdPsi;
/** \brief
* List of all yet created standard zero order assemblers for grdPsi.
*/
static ::std::vector<SubAssembler*> standardSubAssemblersGrdPsi;
static std::vector<SubAssembler*> standardSubAssemblersGrdPsi;
/** \brief
* List of all yet created optimized zero order assemblers for grdPhi.
*/
static ::std::vector<SubAssembler*> optimizedSubAssemblersGrdPhi;
static std::vector<SubAssembler*> optimizedSubAssemblersGrdPhi;
/** \brief
* List of all yet created standard zero order assemblers for grdPhi.
*/
static ::std::vector<SubAssembler*> standardSubAssemblersGrdPhi;
static std::vector<SubAssembler*> standardSubAssemblersGrdPhi;
};
// ============================================================================
......@@ -771,12 +776,12 @@ namespace AMDiS {
/** \brief
* List of all yet created optimized second order assemblers.
*/
static ::std::vector<SubAssembler*> optimizedSubAssemblers;
static std::vector<SubAssembler*> optimizedSubAssemblers;
/** \brief
* List of all yet created standard second order assemblers.
*/
static ::std::vector<SubAssembler*> standardSubAssemblers;
static std::vector<SubAssembler*> standardSubAssemblers;
};
// ============================================================================
......@@ -829,6 +834,8 @@ namespace AMDiS {
*/
Quad2(Operator *op, Assembler *assembler, Quadrature *quad);
~Quad2();
/** \brief
* Implements SubAssembler::calculateElementMatrix().
*/
......@@ -840,6 +847,12 @@ namespace AMDiS {
void calculateElementVector(const ElInfo *, ElementVector */*vec*/) {
ERROR_EXIT("should not be called\n");
};
protected:
/** \brief
* Thread safe temporary vector of DimMats for calculation in calculateElementMatrix().
*/
std::vector< DimMat<double>** > tmpLALt;
};
// ============================================================================
......@@ -886,7 +899,7 @@ namespace AMDiS {
/** \brief
* Thread safe temporary vector of DimMats for calculation in calculateElementMatrix().