small performace speedup for first order assembler.

AMDiS/src/FirstOrderAssembler.cc

@@ -118,25 +118,21 @@ namespace AMDiS {
     Lb.resize(nPoints);
     std::vector<double> phival(nCol);
 
-    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].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 < nCol; i++) {
-	phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));
-      }
+      for (int i = 0; i < nCol; i++)
+	phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));      
 
       for (int i = 0; i < nRow; i++) {
 	(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
-	for (int j = 0; j < nCol; j++) {
-	  mat[i][j] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi) * phival[j];
-	}
+	for (int j = 0; j < nCol; j++)
+	  mat[i][j] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi) * phival[j];	
       }
     }
   }
@@ -150,12 +146,10 @@ namespace AMDiS {
     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].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();
@@ -298,9 +292,8 @@ namespace AMDiS {
 	k = q10->getKVec(i, j);
 	values = q10->getValVec(i, j);
 	double val = 0.0;
-	for (int m = 0; m < nEntries[i][j]; m++) {
-	  val += values[m] * Lb[0][k[m]];
-	}
+	for (int m = 0; m < nEntries[i][j]; m++)
+	  val += values[m] * Lb[0][k[m]];	
 	mat[i][j] += val;
       }
     }
@@ -308,33 +301,33 @@ namespace AMDiS {
 
   Stand01::Stand01(Operator *op, Assembler *assembler, Quadrature *quad) 
     : FirstOrderAssembler(op, assembler, quad, false, GRD_PHI)
-  {}
+  {
+    VectorOfFixVecs<DimVec<double> > 
+      grdPhi(assembler->getRowFESpace()->getMesh()->getDim(), nCol, NO_INIT);
+    tmpGrdPhi.resize(omp_get_overall_max_threads(), grdPhi);
+  }
 
   void Stand01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
   {
-    VectorOfFixVecs<DimVec<double> > grdPhi(dim, nCol, NO_INIT);
-
     const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
     const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
 
     int nPoints = quadrature->getNumPoints();
     int myRank = omp_get_thread_num();
-    VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
+    VectorOfFixVecs<DimVec<double> >& Lb = tmpLb[myRank];
+    VectorOfFixVecs<DimVec<double> >& grdPhi = tmpGrdPhi[myRank];
     Lb.resize(nPoints);
 
-    for (int iq = 0; iq < nPoints; 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++) {
+    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 < nCol; i++) {
+      for (int i = 0; i < nCol; i++)
 	(*(phi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPhi[i]);
-      }
 
       for (int i = 0; i < nRow; i++) {
 	double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
@@ -370,12 +363,10 @@ namespace AMDiS {
     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);
-    }
-    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);
-    }
   
     for (int iq = 0; iq < nPoints; iq++) {
       Lb[iq] *= elInfo->getDet();
@@ -383,10 +374,9 @@ namespace AMDiS {
       const double *psi = psiFast->getPhi(iq);
       grdPhi = phiFast->getGradient(iq);
 
-      for (int i = 0; i < nRow; i++) {
+      for (int i = 0; i < nRow; i++)
 	for (int j = 0; j < nCol; j++)
 	  mat[i][j] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPhi)[j]) * psi[i];
-      }
     }
   }
 

AMDiS/src/FirstOrderAssembler.h

@@ -97,9 +97,13 @@ namespace AMDiS {
     void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat);
 
     /// Implements SubAssembler::calculateElementVector().
-    void calculateElementVector(const ElInfo*, ElementVector&) {
+    void calculateElementVector(const ElInfo*, ElementVector&) 
+    {
       ERROR_EXIT("should not be called\n");
     }
+
+  protected:
+    std::vector<VectorOfFixVecs<DimVec<double> > > tmpGrdPhi;
   };