* Some more optimization stuff

AMDiS/src/Assembler.cc

@@ -1100,32 +1100,41 @@ namespace AMDiS {
 
   Pre2::Pre2(Operator *op, Assembler *assembler, Quadrature *quad) 
     : SecondOrderAssembler(op, assembler, quad, true)
-  {}
+  {
+    q11 = Q11PsiPhi::provideQ11PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
+				      owner->getColFESpace()->getBasisFcts(), 
+				      quadrature);
+    tmpLALt.resize(omp_get_max_threads());
+    for (int i = 0; i < omp_get_max_threads(); i++) {
+      tmpLALt[i] = NEW DimMat<double>*;
+      *(tmpLALt[i]) = NEW DimMat<double>(dim, NO_INIT);
+    }
+  }
+
+  Pre2::~Pre2()
+  {
+    for (int i = 0; i < static_cast<int>(tmpLALt.size()); i++) {
+      DELETE *(tmpLALt[i]);
+      DELETE tmpLALt[i];
+    }
+  }
 
   void Pre2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
   {
-    DimMat<double> **LALt = NEW DimMat<double>*;
-    *LALt=NEW DimMat<double>(dim, NO_INIT);
     const int **nEntries;
     const int *k, *l;
     const double *values;
-    double val;
 
-    if (firstCall) {
-      q11 = Q11PsiPhi::provideQ11PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
-					owner->getColFESpace()->getBasisFcts(), 
-					quadrature);
-      firstCall = false;
-    }
-
-    LALt[0]->set(0.0);
     int myRank = omp_get_thread_num();
+    DimMat<double> **LALt = tmpLALt[myRank];
+    DimMat<double> &tmpMat = *LALt[0];
+    tmpMat.set(0.0);
 
     for (int i = 0; i < static_cast<int>( terms[myRank].size()); i++) {
       (static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, 1, LALt);
     }
 
-    (*LALt[0]) *= elInfo->getDet();
+    tmpMat *= elInfo->getDet();
 
     nEntries = q11->getNumberEntries();
 
@@ -1134,9 +1143,9 @@ namespace AMDiS {
 	k = q11->getKVec(i, i);
 	l = q11->getLVec(i, i);
 	values = q11->getValVec(i, i);
-	val = 0.0;
+	double val = 0.0;
 	for (int m = 0; m < nEntries[i][i]; m++) {
-	  val += values[m] * (*LALt[0])[k[m]][l[m]];
+	  val += values[m] * tmpMat[k[m]][l[m]];
 	}
 	(*mat)[i][i] += val;
 
@@ -1146,7 +1155,7 @@ namespace AMDiS {
 	  values = q11->getValVec(i, j);
 	  val = 0.0;
 	  for (int m = 0; m < nEntries[i][j]; m++) {
-	    val += values[m] * (*LALt[0])[k[m]][l[m]];
+	    val += values[m] * tmpMat[k[m]][l[m]];
 	  }
 	  (*mat)[i][j] += val;
 	  (*mat)[j][i] += val;
@@ -1158,17 +1167,14 @@ namespace AMDiS {
 	  k = q11->getKVec(i, j);
 	  l = q11->getLVec(i, j);
 	  values = q11->getValVec(i, j);
-	  val = 0.0;
+	  double val = 0.0;
 	  for (int m = 0; m < nEntries[i][j]; m++) {
-	    val += values[m] * (*LALt[0])[k[m]][l[m]];
+	    val += values[m] * tmpMat[k[m]][l[m]];
 	  }
 	  (*mat)[i][j] += val;
 	}
       }
     }
-
-    DELETE *LALt;
-    DELETE LALt;
   }
 
   Quad2::Quad2(Operator *op, Assembler *assembler, Quadrature *quad) 

AMDiS/src/Assembler.h

@@ -859,6 +859,11 @@ namespace AMDiS {
      */
     Pre2(Operator *op, Assembler *assembler, Quadrature *quad);
 
+    /** \brief
+     * Destructor.
+     */
+    ~Pre2();
+
     /** \brief
      * Implements SubAssembler::calculateElementMatrix().
      */
@@ -867,7 +872,7 @@ namespace AMDiS {
     /** \brief
      * Implements SubAssembler::calculateElementVector().
      */
-    void calculateElementVector(const ElInfo *, ElementVector */*vec*/) {
+    void calculateElementVector(const ElInfo *, ElementVector *) {
       ERROR_EXIT("should not be called\n");
     };
 
@@ -878,6 +883,11 @@ namespace AMDiS {
      */
     const Q11PsiPhi *q11;
 
+    /** \brief
+     * Thread safe temporary vector of DimMats for calculation in calculateElementMatrix().
+     */
+    ::std::vector< DimMat<double>** > tmpLALt;
+
     friend class SecondOrderAssembler;
   };
 

AMDiS/src/DOFMatrix.cc

@@ -298,33 +298,33 @@ namespace AMDiS {
     if (j!=static_cast<int>(matrix[a].size())) matrix[a][j].col=c;
   }
 
-  double *DOFMatrix::addSparseDOFEntry(double sign, int irow, 
-				       int jcol, double entry,
-				       bool add)
+  void DOFMatrix::addSparseDOFEntry(double sign, int irow, 
+				    int jcol, double entry,
+				    bool add)
   {
     FUNCNAME("DOFMatrix::addSparseDOFEntry()");
 
     MatrixRow *row = &(matrix[irow]);
 
     if (add && !entry) 
-      return NULL;
-
-    double *result = NULL;
+      return;
 
-    int i, freeCol = -1, rowSize = static_cast<int>( row->size());
+    int freeCol = -1;
+    int rowSize = static_cast<int>( row->size());
 
     TEST_EXIT_DBG(jcol >= 0 && 
 		  jcol < colFESpace->getAdmin()->getUsedSize())
       ("Column index %d out of range 0-%d\n", jcol, colFESpace->getAdmin()->getUsedSize() - 1);
 
     // first entry is diagonal entry
-    if (rowFESpace==colFESpace)
+    if (rowFESpace == colFESpace)
       if (rowSize == 0) {
 	MatEntry newEntry = {irow, 0.0};
 	row->push_back(newEntry);
 	rowSize = 1;
       }
 
+    int i;
     // search jcol
     for (i = 0; i < rowSize; i++) {
       // jcol found ?
@@ -350,22 +350,17 @@ namespace AMDiS {
       if (!add) 
 	(*row)[i].entry = 0.0;
       (*row)[i].entry += sign * entry;
-      result = &((*row)[i].entry);
     } else {
-      if(freeCol == -1) {
+      if (freeCol == -1) {
 	MatEntry newEntry = {jcol, sign * entry};
 	row->push_back(newEntry);
-	result = &((*row)[row->size() - 1].entry);
       } else {
 	(*row)[freeCol].col = jcol;
 	if (!add) 
 	  (*row)[freeCol].entry = 0.0;
 	(*row)[freeCol].entry += sign * entry;
-	result = &((*row)[freeCol].entry);
       }
-    }
-  
-    return result;
+    } 
   }
 
   void DOFMatrix::addMatEntry(int row, MatEntry entry)
@@ -527,7 +522,7 @@ namespace AMDiS {
 	      entry += a[rowIndex][j].entry * b[logIndex][physIndex].entry;
 	    }
 	  }
-	  if(entry != 0.0) {
+	  if (entry != 0.0) {
 	    addSparseDOFEntry(1.0, rowIndex, i, entry);
 	  }
 	}
@@ -552,7 +547,7 @@ namespace AMDiS {
 
 	      double entry = aRowIt->entry * bRowIt->entry;
 
-	      if(entry != 0.0) {
+	      if (entry != 0.0) {
 		addSparseDOFEntry(1.0, aCol, bCol, entry);	    
 	      }
 	    }
@@ -592,14 +587,14 @@ namespace AMDiS {
 	// add x contributions to this row
 	for(i=0; i < static_cast<int>((*xIterator).size()); i++) {
 	  colIndex = (*xIterator)[i].col;
-	  if(colIndex >= 0) {
+	  if (colIndex >= 0) {
 	    addSparseDOFEntry(a, rowIndex, colIndex, (*xIterator)[i].entry);
 	  }
 	}
 	// add y contributions to this row
 	for(i=0; i < static_cast<int>((*yIterator).size()); i++) {
 	  colIndex = (*yIterator)[i].col;
-	  if(colIndex >= 0) {
+	  if (colIndex >= 0) {
 	    addSparseDOFEntry(1.0, rowIndex, colIndex, (*yIterator)[i].entry);
 	  }
 	}

AMDiS/src/DOFMatrix.h

@@ -515,9 +515,9 @@ namespace AMDiS {
      * Creates an entry with logical indices irow, icol if there is no entry
      * yet. Than sign * entry is added to the value at this logical indices
      */
-    double *addSparseDOFEntry(double sign,
-			      int irow, int jcol, double entry,
-			      bool add = true);
+    void addSparseDOFEntry(double sign,
+			   int irow, int jcol, double entry,
+			   bool add = true);
 
     inline double *hasSparseDOFEntry(int irow, int jcol) {
       ::std::vector<MatEntry>::iterator it;

AMDiS/src/ElInfo.cc

@@ -178,7 +178,7 @@ namespace AMDiS {
     int dim = mesh_->getDim();
     int posIndex = DIM_OF_INDEX(pos, dim);
     int offset = indexOffset[dim - 1][posIndex];
-  
+    
     return boundary_[offset + i];
   }
 

AMDiS/src/ElInfo.h

@@ -190,7 +190,7 @@ namespace AMDiS {
     /** \brief
      * Get ElInfo's \ref boundary_[i] 
      */
-    virtual BoundaryType getBoundary(int i) const { 
+    inline BoundaryType getBoundary(int i) const { 
       return boundary_[i]; 
     };
 

AMDiS/src/ElInfo3d.cc

@@ -314,8 +314,6 @@ namespace AMDiS {
 
     double det = 0.0;
 
-    int myRank = omp_get_thread_num();
-
     WorldVector<double> *e0 = &tmpWorldVecs[0];
     WorldVector<double> *e1 = &tmpWorldVecs[1];
     WorldVector<double> *e2 = &tmpWorldVecs[2];
@@ -362,13 +360,10 @@ namespace AMDiS {
     const int (*cvg)[4] = NULL;     /* cvg = child_vertex[el_type] */
     int *ce;                    /* ce = child_edge[el_type][ichild] */
     Element *nb, *nbk;
-    const FixVec<Element*, NEIGH> *neigh_old;
     Element *el_old = elinfo_old->element_;
     Flag fillFlag__local = elinfo_old->fillFlag_;
     DegreeOfFreedom *dof;
     int ov = -1;
-    FixVec<Element*, NEIGH> *neigh_local;
-    Flag fill_opp_coords;
     Mesh *mesh = elinfo_old->getMesh();
 
     TEST_EXIT_DBG(el_old->getChild(0))("missing child?\n"); 
@@ -378,7 +373,7 @@ namespace AMDiS {
     fillFlag_ = fillFlag__local;
     parent_  = el_old;
     level_ = elinfo_old->level_ + 1;
-    int el_type_local = ( dynamic_cast<ElInfo3d*>(const_cast<ElInfo*>( elinfo_old)))->getType();
+    int el_type_local = (dynamic_cast<const ElInfo3d*>(elinfo_old))->getType();
     el_type = (el_type_local + 1) % 3;
 
     TEST_EXIT_DBG(element_)("missing child %d?\n", ichild);
@@ -409,8 +404,10 @@ namespace AMDiS {
     if (fillFlag__local.isSet(Mesh::FILL_NEIGH) || 
 	fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
 
-      neigh_local = &neighbour_;
-      neigh_old = &elinfo_old->neighbour_;
+      FixVec<Element*, NEIGH> *neigh_local = &neighbour_;
+      const FixVec<Element*, NEIGH> *neigh_old = &elinfo_old->neighbour_;
+
+      Flag fill_opp_coords;
       fill_opp_coords.setFlags(fillFlag__local & Mesh::FILL_OPP_COORDS);
       
       /*----- nb[0] is other child --------------------------------------------*/

AMDiS/src/Element.h

@@ -110,22 +110,21 @@ namespace AMDiS {
     /** \brief
      * Returns \ref child[0]
      */
-    virtual Element* getFirstChild() const {
+    inline Element* getFirstChild() const {
       return child[0];
     };
 
     /** \brief
      * Returns \ref child[1]
      */
-    virtual Element* getSecondChild() const {
+    inline Element* getSecondChild() const {
       return child[1];
     };
 
     /** \brief
      * Returns \ref child[i], i=0,1
      */
-    virtual Element* getChild(int i) const {
-      FUNCNAME("Element::getChild()");
+    inline Element* getChild(int i) const {
       TEST_EXIT_DBG(i==0 || i==1)("i must be 0 or 1\n");
       return child[i];
     };

AMDiS/src/Lagrange.cc

@@ -790,26 +790,23 @@ namespace AMDiS {
 
     el_info->testFlag(Mesh::FILL_BOUND);
 
-    DimVec<int> parts(dim, NO_INIT);
-    for (int i = 0; i < dim + 1; i++)
-      parts[i] = Global::getGeo(INDEX_OF_DIM(i, dim), dim);
-
-    int index = 0;
-
     // boundaries
+    int index = 0;
     int offset = 0;
     BoundaryType boundaryType;
     for (int i = dim - 1; i > 0; i--)
-      offset += parts[i];
+      offset += Global::getGeo(INDEX_OF_DIM(i, dim), dim);
 
     for (int i = 0; i < dim; i++) {
-      for (int j = offset; j < offset + parts[i]; j++) {
+      int jto = offset + Global::getGeo(INDEX_OF_DIM(i, dim), dim);
+      for (int j = offset; j < jto; j++) {
 	boundaryType = el_info->getBoundary(j);
-	for (int k = 0; k < (*nDOF)[INDEX_OF_DIM(i, dim)]; k++) {
+	int kto = (*nDOF)[INDEX_OF_DIM(i, dim)];
+	for (int k = 0; k < kto; k++) {
 	  result[index++] = boundaryType;
 	}
       }
-      offset -= parts[i+1];
+      offset -= Global::getGeo(INDEX_OF_DIM(i + 1, dim), dim);
     }
 
     // interior nodes in the center

AMDiS/src/OEMSolver.hh

@@ -96,7 +96,7 @@ namespace AMDiS {
     }
 
     if (res <= tolerance) {
-      INFO(info,6)("finished successfully with %d iterations\n", iter);
+	INFO(info,6)("finished successfully with %d iterations\n");
 	return(1);
     }
 

AMDiS/src/Operator.cc

@@ -122,22 +122,23 @@ namespace AMDiS {
   }
 
   void OperatorTerm::l1lt(const DimVec<WorldVector<double> >& Lambda,
-			  DimMat<double>&            LALt,
-			  double                     factor)
+			  DimMat<double>& LALt,
+			  double factor)
   {
-    int                  i, j, k;
-    static const int     dimOfWorld = Global::getGeo(WORLD);
-    int                  dim        = LALt.getNumRows() - 1;
+    static const int dimOfWorld = Global::getGeo(WORLD);
+    int dim = LALt.getNumRows() - 1;
     
     double val = 0.0;
 
-    for (i = 0; i <= dim; i++) {
-      for (val = k = 0; k < dimOfWorld; k++)
+    for (int i = 0; i <= dim; i++) {
+      val = 0.0;
+      for (int k = 0; k < dimOfWorld; k++)
 	val += Lambda[i][k] * Lambda[i][k];
       val *= factor;
       LALt[i][i] += val;
-      for (j = i+1; j <= dim; j++) {
-	for (val = k = 0; k < dimOfWorld; k++)
+      for (int j = i + 1; j <= dim; j++) {
+	val = 0.0;
+	for (int k = 0; k < dimOfWorld; k++)
 	  val += Lambda[i][k] * Lambda[j][k];
 	val *= factor;
 	LALt[i][j] += val;

AMDiS/src/PeriodicBC.cc

@@ -209,9 +209,8 @@ namespace AMDiS {
     TEST_EXIT_DBG(matrix)("no matrix\n");
 
     if (matrix == masterMatrix_) {
-      const BasisFunction *basFcts = rowFESpace->getBasisFcts();
-      int num = basFcts->getNumber();
-      FREE_MEMORY(neighIndices_, DegreeOfFreedom, num);
+      FREE_MEMORY(neighIndices_, DegreeOfFreedom, 
+		  rowFESpace->getBasisFcts()->getNumber());
       masterMatrix_ = NULL;
     }
 
@@ -221,13 +220,15 @@ namespace AMDiS {
     int row, col, newRow, newCol;
     double entry, *newEntryPtr;
 
-    for( rowIt = matrix->begin(), row = 0; rowIt != rowEnd; ++rowIt, ++row) {
+    for (rowIt = matrix->begin(), row = 0; rowIt != rowEnd; ++rowIt, ++row) {
       rowSize = static_cast<int>(rowIt->size());
       newRow = (*associated_)[row];
       for (colIndex = 0; colIndex < rowSize; colIndex++) {
 	col = (*rowIt)[colIndex].col;
-	if(col == DOFMatrix::UNUSED_ENTRY) continue;
-	if(col == DOFMatrix::NO_MORE_ENTRIES) break;
+	if (col == DOFMatrix::UNUSED_ENTRY) 
+	  continue;
+	if (col == DOFMatrix::NO_MORE_ENTRIES) 
+	  break;
 	newCol = (*associated_)[col];
 
 	newEntryPtr = matrix->hasSparseDOFEntry(newRow, newCol);

AMDiS/src/Preconditioner.h

@@ -241,7 +241,7 @@ namespace AMDiS {
       int size = scalPrecons.getSize();
 
 #ifdef _OPENMP
-#pragma omp parallel for
+#pragma omp parallel for num_threads(size)
 #endif
       for (i = 0; i < size; i++) {
 	scalPrecons[i]->init();
@@ -255,7 +255,7 @@ namespace AMDiS {
       int i;
       int size = scalPrecons.getSize();
 #ifdef _OPENMP
-#pragma omp parallel for
+#pragma omp parallel for num_threads(size)
 #endif
       for (i = 0; i < size; i++) {
 	scalPrecons[i]->precon(x->getDOFVector(i));
@@ -269,7 +269,7 @@ namespace AMDiS {
       int i;
       int size = scalPrecons.getSize();
 #ifdef _OPENMP
-#pragma omp parallel for
+#pragma omp parallel for num_threads(size)
 #endif
       for (i = 0; i < size; i++) {
 	scalPrecons[i]->exit();

AMDiS/src/ProblemVec.cc

@@ -701,6 +701,8 @@ namespace AMDiS {
 #pragma omp parallel for 
 #endif
     for (i = 0; i < numComponents_; i++) {
+      const BasisFunction *basisFcts = componentSpaces_[i]->getBasisFcts();
+
       for (int j = 0; j < numComponents_; j++) {
 	// Only if this variable is true, the current matrix will be assembled.	
 	bool assembleMatrix = true;
@@ -734,7 +736,7 @@ namespace AMDiS {
 	} else {
 	  BoundaryType *bound = NULL;
 	  if (useGetBound_) {
-	    bound = GET_MEMORY(BoundaryType, componentSpaces_[i]->getBasisFcts()->getNumber());
+	    bound = GET_MEMORY(BoundaryType, basisFcts->getNumber());
 	  }
 	  
 	  TraverseStack stack;
@@ -742,7 +744,7 @@ namespace AMDiS {
 	  
 	  while (elInfo) {
 	    if (useGetBound_) {
-	      componentSpaces_[i]->getBasisFcts()->getBound(elInfo, bound);
+	      basisFcts->getBound(elInfo, bound);
 	    }
 	          
 	    if (assembleMatrix) {
@@ -767,7 +769,7 @@ namespace AMDiS {
 	  	    matrix->getBoundaryManager()->exitMatrix(matrix);	  
 
 	  if (useGetBound_) {
-	    FREE_MEMORY(bound, BoundaryType, componentSpaces_[i]->getBasisFcts()->getNumber());
+	    FREE_MEMORY(bound, BoundaryType, basisFcts->getNumber());
 	  }	  
 	}
 	

AMDiS/src/RefinementManager.cc

@@ -15,9 +15,7 @@ namespace AMDiS {
   RefinementManager* RefinementManager::traversePtr = NULL;
   bool RefinementManager::doMoreRecursiveRefine = false;
   int RefinementManager::callRefineInterpol = 0;
-
   RCNeighbourList* RefinementManager3d::refList = NULL;
-
   int RefinementManager::globalMark = 0;
 
   int RefinementManager::globalRefineFct(ElInfo* elinfo)
@@ -32,7 +30,9 @@ namespace AMDiS {
   {
     FUNCNAME("RefinementManager::globalRefine()");
 
-    if(mark <= 0) return static_cast<Flag>(0);
+    if (mark <= 0) 
+      return static_cast<Flag>(0);
+
     globalMark = mark;
     aMesh->traverse(-1, 
 		    Mesh::CALL_LEAF_EL |
@@ -50,7 +50,6 @@ namespace AMDiS {
     FUNCNAME("RefinementManager::refineMesh()");
 
     mesh = aMesh;
-
     int n_elements = mesh->getNumberOfLeaves();
     ElInfo *el_info;
 

AMDiS/src/RefinementManager2d.cc

@@ -17,8 +17,9 @@ namespace AMDiS {
 
   ElInfo* RefinementManager2d::refineFunction(ElInfo* el_info)
   {
-    FUNCNAME("RefinementManager::refineFunction");
-    int  n_neigh;
+    FUNCNAME("RefinementManager::refineFunction()");
+
+    int n_neigh;
     bool bound = false;
 
     DegreeOfFreedom *edge[2];
@@ -76,7 +77,7 @@ namespace AMDiS {
     ::std::map<int, VertexVector*>::iterator it;
     ::std::map<int, VertexVector*>::iterator end = mesh->getPeriodicAssociations().end();
 
-    while(edge[0] != NULL) {
+    while (edge[0] != NULL) {
       periodicList = refineList->periodicSplit(edge, 
 					       next_edge,
 					       &n_neigh,
@@ -134,8 +135,6 @@ namespace AMDiS {
     /*  and now refine the patch                                                */
     /****************************************************************************/
 
-    //newDOF = refinePatch(edge, refineList, n_neigh, bound);
-
     DELETE refineList;
 
     return(el_info);
@@ -206,63 +205,10 @@ namespace AMDiS {
     /****************************************************************************/
     /*  first refine the element                                                */
     /****************************************************************************/
-  
-    // check for periodic boundary
-    //   DegreeOfFreedom *periodicDOF[3] = {NULL, NULL, NULL};
-    //   RCNeighbourList *periodicRefineList = NULL;
-    //   int n_neigh_periodic;
-
-    //   if (neigh && 
-    //       (neigh->getDOF(0) != el->getDOF(1)) && 
-    //       (neigh->getDOF(0) != el->getDOF(0))) 
-    //   {
-    //     periodicDOF[0] = mesh->getDOF(VERTEX);
-    //     mesh->incrementNumberOfVertices(1);
-    //     mesh->incrementNumberOfEdges(1);
-    //     if (mesh->getNumberOfDOFs(EDGE)) {
-    //       periodicDOF[1] = mesh->getDOF(EDGE);
-    //       periodicDOF[2] = mesh->getDOF(EDGE);
-    //     }
-
-    //     DegreeOfFreedom *edge[2] = { 
-    //       const_cast<DegreeOfFreedom*>(el->getDOF(0)), 
-    //       const_cast<DegreeOfFreedom*>(el->getDOF(1)) 
-    //     };
-    //     DegreeOfFreedom *periodicEdge[2] = { 
-    //       const_cast<DegreeOfFreedom*>(neigh->getDOF(0)), 
-    //       const_cast<DegreeOfFreedom*>(neigh->getDOF(1)) 
-    //     };
-
-    //     periodicRefineList = refineList->periodicSplit(edge,
-    // 						   &n_neigh,
-    // 						   periodicEdge,
-    // 						   &n_neigh_periodic);
-
-    // //     dynamic_cast<LeafDataPeriodic*>(neigh->getElementData(PERIODIC))->
-    // //       setNewDOF(dof[0][0]);
-    // //     dynamic_cast<LeafDataPeriodic*>(el->getElementData(PERIODIC))->
-    // //       setNewDOF(periodicDOF[0][0]);
-
-    // //     MSG("el %d new dof %d periodic dof %d\n",
-    // // 	el->getIndex(),
-    // // 	dof[0][0],
-    // // 	periodicDOF[0][0]);