* Added "store symbolic" parameter for UMFPACK-solver to speed up solving...

* Added "store symbolic" parameter for UMFPACK-solver to speed up solving systems resulting from non-adaptive meshes

AMDiS/src/MacroReader.cc

@@ -475,7 +475,7 @@ namespace AMDiS {
     FUNCNAME("MacroInfo::readAMDiSMacro");
     FILE       *file;
     int        dim;
-    int        dow, nv, ne, i, j, k;
+    int        dow, nv, ne, j, k;
     double     dbl;
     char       name[128], line[256];
     int        line_no, n_keys, i_key, sort_key[N_KEYS], nv_key, ne_key;
@@ -600,15 +600,13 @@ namespace AMDiS {
 	    break;
 	  case 4:
 	    fscanf(file, "%*s %*s");
-	    for (i = 0; i < nv; i++) 
-	      {
-		for (j = 0; j <Global::getGeo(WORLD) ; j++)
-		  {
-		    TEST_EXIT(fscanf(file, "%lf", &dbl) == 1)
-		      ("error while reading coordinates, check file %s\n", name);
-		    coords[i][j] = dbl;
-		  }
+	    for (int i = 0; i < nv; i++) {
+	      for (j = 0; j <Global::getGeo(WORLD) ; j++) {
+		TEST_EXIT(fscanf(file, "%lf", &dbl) == 1)
+		  ("error while reading coordinates, check file %s\n", name);
+		coords[i][j] = dbl;
 	      }
+	    }
 	    key_def[4] = true;
 	    break;
 	  case 5:
@@ -617,7 +615,7 @@ namespace AMDiS {
 	    /* global number of vertices on a single element                            */
 	    /****************************************************************************/
 
-	    for (i = 0; i < ne; i++) {
+	    for (int i = 0; i < ne; i++) {
 		TEST_EXIT(read_indices(file, *ind))
 		  ("cannot read vertex indices of element %d in file %s\n",  i, name);
 
@@ -632,8 +630,7 @@ namespace AMDiS {
 	    /****************************************************************************/
 	    /* MEL boundary pointers                                                    */
 	    /****************************************************************************/
-
-	    for (i = 0; i < ne; i++) {  
+	    for (int i = 0; i < ne; i++) {  
                 // boundary information of ith element 
 
 		TEST_EXIT(read_indices(file, *ind))
@@ -658,7 +655,7 @@ namespace AMDiS {
 	    /*   else init them by a call of fill_mel_neighbour()                       */
 	    /****************************************************************************/
 	    neigh_set = true;
-	    for (i = 0; i < ne; i++) {
+	    for (int i = 0; i < ne; i++) {
                 //  neighbour information about ith element
 
 		if (read_indices(file, *ind))
@@ -682,7 +679,7 @@ namespace AMDiS {
 	    if (dim == 2 || dim == 1)
 	      ERROR("there is no element type in 2d and 2d; ignoring data for elType\n");
 
-	    for (i = 0; i < ne; i++) {
+	    for (int i = 0; i < ne; i++) {
 		TEST_EXIT(fscanf(file, "%d", &j) == 1)
 		  ("cannot read elType of element %d in file %s\n",
 		   i, name);
@@ -700,11 +697,11 @@ namespace AMDiS {
 	      int numFaces = mesh->getGeo(FACE);
 	      int numEdgesAtBoundary = 0;
 
-	      for(k = 1; k < dim; k++) {
+	      for (k = 1; k < dim; k++) {
 		numEdgesAtBoundary += k;
 	      }
 
-	      for (i = 0; i < ne; i++) {
+	      for (int i = 0; i < ne; i++) {
 		TEST_EXIT(read_indices(file, *ind))
 		  ("cannot read boundary projector of element %d in file %s\n", i, name);
 	
@@ -736,7 +733,7 @@ namespace AMDiS {
 	    /* MEL regions                                                              */
 	    /****************************************************************************/
 
-	    for (i = 0; i < ne; i++) {
+	    for (int i = 0; i < ne; i++) {
 	      TEST_EXIT(fscanf(file, "%d", &j) == 1)
 		("cannot read region of element %d in file %s\n", i, name);
 	      if(j >= 0) {
@@ -751,7 +748,7 @@ namespace AMDiS {
 	    break;
 	  case 11:
 	    fscanf(file, "%*s %*s");
-	    for (i = 0; i < ne; i++) {
+	    for (int i = 0; i < ne; i++) {
 	      TEST_EXIT(read_indices(file, *ind))
 		("cannot read surface regions of element %d in file %s\n", i, name);
 
@@ -807,10 +804,8 @@ namespace AMDiS {
 
   void MacroInfo::dirichletBoundary()
   {
-    int      i, k;
-
-    for (i = 0; i < static_cast<int>( mel.size()); i++) {
-      for (k = 0; k < mesh->getGeo(NEIGH); k++) {
+    for (int i = 0; i < static_cast<int>( mel.size()); i++) {
+      for (int k = 0; k < mesh->getGeo(NEIGH); k++) {
 	if (mel[i]->neighbour[k])
 	  mel[i]->boundary[k] = INTERIOR;
 	else

AMDiS/src/OEMSolver.h

@@ -135,7 +135,9 @@ namespace AMDiS {
     /** \brief
      * Returns solvers \ref name.
      */
-    inline const ::std::string& getName() { return name; };
+    inline const ::std::string& getName() { 
+      return name; 
+    };
 
     /** \brief
      * Returns \ref tolerance
@@ -158,10 +160,6 @@ namespace AMDiS {
       return residual;
     };
 
-    //   inline Preconditioner<VectorType> *getLeftPrecon() { return leftPrecon; };
-
-    //   inline Preconditioner<VectorType> *getRightPrecon() { return rightPrecon; };
-
     /** \} */ 
 
     // ===== setting-methods ======================================================
@@ -173,32 +171,28 @@ namespace AMDiS {
     /** \brief
      * Sets \ref tolerance
      */
-    inline void setTolerance(double tol) 
-    {
+    inline void setTolerance(double tol) {
       tolerance = tol;
     };
 
     /** \brief
      * Sets \ref relative
      */
-    inline void setRelative(bool rel)
-    {
+    inline void setRelative(bool rel) {
       relative = rel;
     };
 
     /** \brief
      * Sets \ref max_iter
      */
-    inline void setMaxIterations(int i) 
-    {
+    inline void setMaxIterations(int i) {
       max_iter = i;
     };
 
     /** \brief
      * Sets \ref info
      */
-    inline void setInfo(int i)
-    {
+    inline void setInfo(int i) {
       info = i;
     };
   
@@ -211,14 +205,6 @@ namespace AMDiS {
       return vectorCreator;
     };
 
-    //   inline void setLeftPrecon(Preconditioner<VectorType> *precon) {
-    //     leftPrecon = precon;
-    //   };
-
-    //   inline void setRightPrecon(Preconditioner<VectorType> *precon) {
-    //     rightPrecon = precon;
-    //   };
-
     /** \} */
 
   protected:

AMDiS/src/ProblemVec.cc

@@ -192,19 +192,19 @@ namespace AMDiS {
     GET_PARAMETER(0, name_ + "->mesh", &meshName);
     TEST_EXIT(meshName != "")("no mesh name spezified\n");
     int dim = 0;
+    ::std::cout << name_ + "->dim" << ::std::endl;
     GET_PARAMETER(0, name_ + "->dim", "%d", &dim);
     TEST_EXIT(dim)("no problem dimension spezified!\n");
 
-    int i, refSet;
-    for(i = 0; i < numComponents_; i++) {
+    for (int i = 0; i < numComponents_; i++) {
       sprintf(number, "%d", i);
-      refSet = -1;
-      GET_PARAMETER(0, name_ + "->refinement set[" + number + "]","%d", &refSet);
-      if(refSet < 0) {
+      int refSet = -1;
+      GET_PARAMETER(0, name_ + "->refinement set[" + number + "]", "%d", &refSet);
+      if (refSet < 0) {
 	WARNING("refinement set for component %d not set\n", i);
 	refSet = 0;
       }
-      if(meshForRefinementSet[refSet] == NULL) {
+      if (meshForRefinementSet[refSet] == NULL) {
 	Mesh *newMesh = NEW Mesh(meshName, dim);
 	meshForRefinementSet[refSet] = newMesh;
 	meshes_.push_back(newMesh);
@@ -237,13 +237,13 @@ namespace AMDiS {
     char number[3];
 
     ::std::map< ::std::pair<Mesh*, int>, FiniteElemSpace*> feSpaceMap;
-    int i, dim = -1;
+    int dim = -1;
     GET_PARAMETER(0, name_ + "->dim", "%d", &dim);
     TEST_EXIT(dim != -1)("no problem dimension spezified!\n");
 
     componentSpaces_.resize(numComponents_, NULL);
 
-    for(i = 0; i < numComponents_; i++) {
+    for (int i = 0; i < numComponents_; i++) {
       sprintf(number, "%d", i);
       GET_PARAMETER(0, name_ + "->polynomial degree[" + number + "]","%d", &degree);
 
@@ -263,7 +263,7 @@ namespace AMDiS {
     }
 
     // create dof admin for vertex dofs if neccessary
-    for (i = 0; i < static_cast<int>(meshes_.size()); i++) {
+    for (int i = 0; i < static_cast<int>(meshes_.size()); i++) {
       if (meshes_[i]->getNumberOfDOFs(VERTEX) == 0) {
 	DimVec<int> ln_dof(meshes_[i]->getDim(), DEFAULT_VALUE, 0);
 	ln_dof[VERTEX]= 1;
@@ -892,13 +892,10 @@ namespace AMDiS {
 
     PeriodicBC *periodic = new PeriodicBC(type, feSpace);
 
-    //   TEST_EXIT(mesh->getPeriodicBCMap()[type] == NULL)
-    //     ("periodic condition already set in mesh\n");
-    //   mesh->getPeriodicBCMap()[type] = periodic;
-
-    if(systemMatrix_ && (*systemMatrix_)[row][col]) 
+    if (systemMatrix_ && (*systemMatrix_)[row][col]) 
       (*systemMatrix_)[row][col]->getBoundaryManager()->addBoundaryCondition(periodic);
-    if(rhs_) 
+
+    if (rhs_) 
       rhs_->getDOFVector(row)->getBoundaryManager()->
 	addBoundaryCondition(periodic);
   }

AMDiS/src/UmfPackSolver.h

@@ -102,6 +102,20 @@ namespace AMDiS {
      * These vectors are justed to calculate the final residual of the solution.
      */
     VectorType *r, *p;
+
+    /** \brief
+     * Stores the result of umfpack_di_symbolic (the symbolic analysis which is needed
+     * for the numeric factorization).
+     */
+    void *symbolic_;
+
+    /** \brief
+     * If the symbolic analysis should be done only once (for example, if the matrices
+     * to solve have all the same pattern because of no adaptivity), this variable 
+     * is one, otherwise 0 and the symbolic analysis will be performed at every call
+     * of solveSystem().
+     */
+    int storeSymbolic_;
   };
   
 }

AMDiS/src/UmfPackSolver.hh

@@ -11,11 +11,20 @@ namespace AMDiS {
 
   template<typename VectorType>
   UmfPackSolver<VectorType>::UmfPackSolver(::std::string name) 
-    : OEMSolver<VectorType>(name) 
-  {}
+    : OEMSolver<VectorType>(name),
+      symbolic_(NULL),
+      storeSymbolic_(0)
+  {
+    GET_PARAMETER(0, name + "->store symbolic", "%d", &storeSymbolic_);
+  }
 
   template<typename VectorType>
-  UmfPackSolver<VectorType>::~UmfPackSolver() {}
+  UmfPackSolver<VectorType>::~UmfPackSolver() 
+  {
+    if (storeSymbolic_ && (symbolic_ != NULL)) {
+      umfpack_di_free_symbolic(&symbolic_);
+    }
+  }
 
   template<typename VectorType>
   int UmfPackSolver<VectorType>::solveSystem(MatVecMultiplier<VectorType> *matVec,
@@ -105,7 +114,7 @@ namespace AMDiS {
     }
 
 	   
-    void *Symbolic, *Numeric;
+    void *numeric;
     double Control[UMFPACK_CONTROL];
     double Info[UMFPACK_INFO];
     int status = 0;
@@ -116,24 +125,32 @@ namespace AMDiS {
     umfpack_di_defaults(Control);
    
     // Run UMFPack.
-    status = umfpack_di_symbolic(newMatrixSize, newMatrixSize, Ap, Ai, Ax, &Symbolic, Control, Info);   
-    if (!status == UMFPACK_OK) {
-      ERROR_EXIT("UMFPACK Error in function umfpack_di_symbolic");
-    }    
- 
-    status = umfpack_di_numeric(Ap, Ai, Ax, Symbolic, &Numeric, Control, Info);
-    if (!status == UMFPACK_OK) {
-      ERROR_EXIT("UMFPACK Error in function umfpack_di_numberic");
+    if (!(storeSymbolic_) || 
+	(storeSymbolic_ && (symbolic_ == NULL))) {
+
+      status = umfpack_di_symbolic(newMatrixSize, newMatrixSize, Ap, Ai, Ax, &symbolic_, Control, Info);   
+      if (!status == UMFPACK_OK) {
+	ERROR_EXIT("UMFPACK Error in function umfpack_di_symbolic");
+      }  
+      ::std::cout << "symbolic durchgefuehrt" << ::std::endl;
+    } else {
+      ::std::cout << "symbolic ausgelassen" << ::std::endl;
     }
     
-    status = umfpack_di_solve(UMFPACK_A, Ap, Ai, Ax, xvec, bvec, Numeric, Control, Info);
+    status = umfpack_di_numeric(Ap, Ai, Ax, symbolic_, &numeric, Control, Info);
+    if (!status == UMFPACK_OK) {
+      ERROR_EXIT("UMFPACK Error in function umfpack_di_numeric");
+    }
+  
+    status = umfpack_di_solve(UMFPACK_A, Ap, Ai, Ax, xvec, bvec, numeric, Control, Info);
     if (!status == UMFPACK_OK) {
       ERROR_EXIT("UMFPACK Error in function umfpack_di_solve");
     }
-
-    
-    umfpack_di_free_symbolic(&Symbolic);
-    umfpack_di_free_numeric(&Numeric);
+   
+    if (!storeSymbolic_) {
+      umfpack_di_free_symbolic(&symbolic_);
+    }
+    umfpack_di_free_numeric(&numeric);
 
     
     // Copy and resort solution.