Fixed nnz caluculations for petsc solver.

AMDiS/src/parallel/InteriorBoundary.cc

@@ -228,7 +228,14 @@ namespace AMDiS {
 		      bound.rankObj.elIndex)("Should not happen!\n");
 
 	bound.rankObj.el = elIndexMap[bound.rankObj.elIndex];
-	bound.neighObj.el = NULL;
+
+	// For the case of periodic interior boundaries, a rank may have an boundary
+	// with itself. In this case, also the pointer to the neighbour object must
+	// be set correctly.
+	if (elIndexMap.count(bound.neighObj.elIndex))
+	  bound.neighObj.el = elIndexMap[bound.neighObj.elIndex];
+	else
+	  bound.neighObj.el = NULL;
       }
     }
   }

AMDiS/src/parallel/MeshDistributor.cc

@@ -1644,6 +1644,10 @@ namespace AMDiS {
 
 	for (unsigned int i = 0; i < it->second.size(); i++) {
 	  AtomicBoundary &bound = it->second[i];
+
+	  TEST_EXIT_DBG(bound.rankObj.el)("No rank object!\n");
+	  TEST_EXIT_DBG(bound.neighObj.el)("No neigh object!\n");
+
 	  DofContainer dofs0, dofs1;
 	  bound.rankObj.el->getVertexDofs(feSpace, bound.rankObj, dofs0);
 	  bound.rankObj.el->getNonVertexDofs(feSpace, bound.rankObj, dofs0); 

AMDiS/src/parallel/MeshDistributor.h

@@ -201,7 +201,10 @@ namespace AMDiS {
     /// rank's partition, but is owned by some other rank.
     inline bool getIsRankDof(DegreeOfFreedom dof)
     {
-      return isRankDof[dof];
+      if (isRankDof.count(dof))
+	return isRankDof[dof];
+
+      return false;
     }
 
     inline long getLastMeshChangeIndex()
@@ -526,7 +529,7 @@ namespace AMDiS {
      * If periodic boundaries are used, this map stores to each periodic DOF in rank's
      * partition the set of periodic boundaries the DOF is associated to. In 2D, most
      * DOFs are only on one periodic boundary. Only, e.g., in a box with all boundaries
-     * being periodic, the for corners are associated by two different boundaries.     
+     * being periodic, the four corners are associated by two different boundaries.
      */
     map<int, std::set<BoundaryType> > periodicDofAssociations;
 

AMDiS/src/parallel/ParallelDebug.cc

@@ -143,32 +143,6 @@ namespace AMDiS {
   {
     FUNCNAME("ParallelDebug::testPeriodicBoundary()");
 
-    DegreeOfFreedom testDof = pdb.mapGlobalToLocal(968);
-    if (testDof != -1) {
-      WorldVector<double> c;
-      pdb.mesh->getDofIndexCoords(testDof, pdb.feSpace, c);
-
-      MSG("FOUND: LOCAL %d   GLOBAL %d  COORDS %f %f %f\n",
-	  testDof, 968, c[0], c[1], c[2]);
-    }
-
-    testDof = pdb.mapGlobalToLocal(973);
-    if (testDof != -1) {
-      WorldVector<double> c;
-      pdb.mesh->getDofIndexCoords(testDof, pdb.feSpace, c);
-
-      MSG("FOUND: LOCAL %d   GLOBAL %d  COORDS %f %f %f\n",
-	  testDof, 973, c[0], c[1], c[2]);
-    }
-
-    testDof = pdb.mapGlobalToLocal(1795);
-    if (testDof != -1) {
-      WorldVector<double> c;
-      pdb.mesh->getDofIndexCoords(testDof, pdb.feSpace, c);
-
-      MSG("FOUND: LOCAL %d   GLOBAL %d  COORDS %f %f %f\n",
-	  testDof, 1795, c[0], c[1], c[2]);
-    }
     // === 1. check: All periodic DOFs should have at least a correct number ===
     // === of periodic associations.                                         ===
 

AMDiS/src/parallel/PetscSolver.cc

@@ -100,6 +100,9 @@ namespace AMDiS {
 	  // Test if the current col dof is a periodic dof.
 	  bool periodicCol = meshDistributor->isPeriodicDof(globalColDof);
 
+	  if (value(*icursor) == 0.0 && globalRowDof != globalColDof)
+	    continue;
+
 	  if (!periodicCol) {
 	    // Calculate the exact position of the column index in the petsc matrix.
 	    cols.push_back(globalColDof * dispMult + dispAddCol);
@@ -130,7 +133,7 @@ namespace AMDiS {
 	      }
 	    }
 
-	    for (int i = 0; i < newCols.size(); i++) {
+	    for (unsigned int i = 0; i < newCols.size(); i++) {
 	      cols.push_back(newCols[i] * dispMult + dispAddCol);
 	      values.push_back(scaledValue);	      
 	    }
@@ -140,12 +143,18 @@ namespace AMDiS {
 	MatSetValues(petscMatrix, 1, &rowIndex, cols.size(), 
 		     &(cols[0]), &(values[0]), ADD_VALUES);	
       } else {
+	std::map<int, std::vector<int> > colsMap;
+	std::map<int, std::vector<double> > valsMap;
+
 	for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
-	     icursor != icend; ++icursor) {	  
+	     icursor != icend; ++icursor) {
 
 	  // Global index of the current column index.
 	  int globalColDof = meshDistributor->mapLocalToGlobal(col(*icursor));
 
+	  if (value(*icursor) == 0.0 && globalRowDof != globalColDof)
+	    continue;
+
 	  std::set<int> perAsc;
 
 	  if (meshDistributor->isPeriodicDof(globalColDof)) {
@@ -157,14 +166,12 @@ namespace AMDiS {
 		perAsc.insert(*it);
 	  }
 
-	  if (meshDistributor->isPeriodicDof(globalRowDof)) {
-	    std::set<int>& perRowAsc = 
-	      meshDistributor->getPerDofAssociations(globalRowDof);
-	    for (std::set<int>::iterator it = perRowAsc.begin(); 
-		 it != perRowAsc.end(); ++it)
-	      if (*it >= -3)
-		perAsc.insert(*it);
-	  }
+	  std::set<int>& perRowAsc = 
+	    meshDistributor->getPerDofAssociations(globalRowDof);
+	  for (std::set<int>::iterator it = perRowAsc.begin(); 
+	       it != perRowAsc.end(); ++it)
+	    if (*it >= -3)
+	      perAsc.insert(*it);
 
 	  double scaledValue = 
 	    value(*icursor) * std::pow(0.5, static_cast<double>(perAsc.size()));
@@ -195,10 +202,23 @@ namespace AMDiS {
 	       eIt != entry.end(); ++eIt) {
 	    // Calculate petsc row index.
 	    int rowIndex = eIt->first * dispMult + dispAddRow;
+
 	    int colIndex = eIt->second * dispMult + dispAddCol;
-	    MatSetValue(petscMatrix, rowIndex, colIndex, scaledValue, ADD_VALUES);
+	    colsMap[rowIndex].push_back(colIndex);
+	    valsMap[rowIndex].push_back(scaledValue);
 	  }
 	}
+
+
+	for (std::map<int, std::vector<int> >::iterator rowIt = colsMap.begin();
+	     rowIt != colsMap.end(); ++rowIt) {
+	  TEST_EXIT_DBG(rowIt->second.size() == valsMap[rowIt->first].size())
+	    ("Should not happen!\n");
+
+	  int rowIndex = rowIt->first;
+	  MatSetValues(petscMatrix, 1, &rowIndex, rowIt->second.size(),
+		       &(rowIt->second[0]), &(valsMap[rowIt->first][0]), ADD_VALUES);
+	}
       }
     }
   }
@@ -256,111 +276,123 @@ namespace AMDiS {
     namespace traits = mtl::traits;
     typedef DOFMatrix::base_matrix_type Matrix;
     typedef std::vector<std::pair<int, int> > MatrixNnzEntry;
+    typedef std::map<int, DofContainer> RankToDofContainer;
 
     // Stores to each rank a list of nnz entries (i.e. pairs of row and column index)
-    // that this rank will send to. This nnz entries will be assembled on this rank,
+    // that this rank will send to. These nnz entries will be assembled on this rank,
     // but because the row DOFs are not DOFs of this rank they will be send to the
     // owner of the row DOFs.
     std::map<int, MatrixNnzEntry> sendMatrixEntry;
 
+    // Maps to each DOF that must be send to another rank the rank number of the
+    // receiving rank.
+    std::map<DegreeOfFreedom, int> sendDofToRank;
+
 
     // First, create for all ranks we send data to MatrixNnzEntry object with 0 entries.
-    typedef std::map<int, DofContainer> RankToDofContainer;
-    RankToDofContainer& recvDofs = meshDistributor->getRecvDofs();
-    for (RankToDofContainer::iterator it = recvDofs.begin();
-	 it != recvDofs.end(); ++it)
+    for (RankToDofContainer::iterator it = meshDistributor->getRecvDofs().begin();
+	 it != meshDistributor->getRecvDofs().end(); ++it) {
       sendMatrixEntry[it->first].resize(0);
 
+      for (DofContainer::iterator dofIt = it->second.begin();
+	   dofIt != it->second.end(); ++dofIt)
+	sendDofToRank[**dofIt] = it->first;
+    }
+
+
+    std::set<int> recvFromRank;
+    for (RankToDofContainer::iterator it = meshDistributor->getSendDofs().begin();
+	 it != meshDistributor->getSendDofs().end(); ++it)
+      recvFromRank.insert(it->first);
+
 
     for (int i = 0; i < nComponents; i++) {
       for (int j = 0; j < nComponents; j++) {
- 	if ((*mat)[i][j]) {
-	  Matrix bmat = (*mat)[i][j]->getBaseMatrix();
+ 	if (!(*mat)[i][j])
+	  continue;
+
+	Matrix bmat = (*mat)[i][j]->getBaseMatrix();
 
-	  traits::col<Matrix>::type col(bmat);
-	  traits::const_value<Matrix>::type value(bmat);
+	traits::col<Matrix>::type col(bmat);
+	traits::const_value<Matrix>::type value(bmat);
 	  
-	  typedef traits::range_generator<row, Matrix>::type cursor_type;
-	  typedef traits::range_generator<nz, cursor_type>::type icursor_type;
+	typedef traits::range_generator<row, Matrix>::type cursor_type;
+	typedef traits::range_generator<nz, cursor_type>::type icursor_type;
+	
+	for (cursor_type cursor = begin<row>(bmat), 
+	       cend = end<row>(bmat); cursor != cend; ++cursor) {
 	  
-	  for (cursor_type cursor = begin<row>(bmat), 
-		 cend = end<row>(bmat); cursor != cend; ++cursor) {
-
-	    // Map the local row number to the global DOF index and create from it
-	    // the global PETSc row index of this DOF.
-	    int petscRowIdx = 
-	      meshDistributor->mapLocalToGlobal(*cursor) * nComponents + i;
-
-	    if (meshDistributor->getIsRankDof(*cursor)) {
+	  int globalRowDof = meshDistributor->mapLocalToGlobal(*cursor);
 
-	      // === The current row DOF is a rank dof, so create the corresponding ===
-	      // === nnz values directly on rank's nnz data.                        ===
-
-	      // This is the local row index of the local PETSc matrix.
-	      int localPetscRowIdx = 
-		petscRowIdx - meshDistributor->getRstart() * nComponents;
+	  std::vector<int> rows;
+	  rows.push_back(globalRowDof);
+	  std::vector<int> rowRank;
+	  if (meshDistributor->getIsRankDof(*cursor)) {
+	    rowRank.push_back(meshDistributor->getMpiRank());
+	  } else {
+	    // Find out who is the member of this DOF.
+	    TEST_EXIT_DBG(sendDofToRank.count(*cursor))("Should not happen!\n");
+	    
+	    rowRank.push_back(sendDofToRank[*cursor]);
+	  }
 
-	      TEST_EXIT_DBG(localPetscRowIdx >= 0 && localPetscRowIdx < nRankRows)
-		("Should not happen! \n Debug info: localRowIdx = %d   globalRowIndx = %d   petscRowIdx = %d   localPetscRowIdx = %d   rStart = %d   nCompontens = %d   nRankRows = %d\n",
-		 *cursor, meshDistributor->mapLocalToGlobal(*cursor), petscRowIdx, localPetscRowIdx, meshDistributor->getRstart(), nComponents, nRankRows);
+	  // Map the local row number to the global DOF index and create from it
+	  // the global PETSc row index of this DOF.
+	  
+	  int petscRowIdx = globalRowDof * nComponents + i;
+	  
+	  if (meshDistributor->getIsRankDof(*cursor)) {
+	    	    
+	    // === The current row DOF is a rank dof, so create the corresponding ===
+	    // === nnz values directly on rank's nnz data.                        ===
+	    
+	    // This is the local row index of the local PETSc matrix.
+	    int localPetscRowIdx = 
+	      petscRowIdx - meshDistributor->getRstart() * nComponents;
+	    
+	    TEST_EXIT_DBG(localPetscRowIdx >= 0 && localPetscRowIdx < nRankRows)
+	      ("Should not happen! \n Debug info: localRowIdx = %d   globalRowIndx = %d   petscRowIdx = %d   localPetscRowIdx = %d   rStart = %d   nCompontens = %d   nRankRows = %d\n",
+	       *cursor, meshDistributor->mapLocalToGlobal(*cursor), petscRowIdx, localPetscRowIdx, meshDistributor->getRstart(), nComponents, nRankRows);
+	    
+	    
+	    // Traverse all non zero entries in this row.
+	    for (icursor_type icursor = begin<nz>(cursor), 
+		   icend = end<nz>(cursor); icursor != icend; ++icursor) {
+	      int petscColIdx = 
+		meshDistributor->mapLocalToGlobal(col(*icursor)) * nComponents + j;
 	      
-	      // Traverse all non zero entries in this row.
-	      for (icursor_type icursor = begin<nz>(cursor), 
-		     icend = end<nz>(cursor); icursor != icend; ++icursor) {
+	      if (value(*icursor) != 0.0 || petscRowIdx == petscColIdx) {
+		// The row DOF is a rank DOF, if also the column is a rank DOF, 
+		// increment the d_nnz values for this row, otherwise the o_nnz value.
+		if (petscColIdx >= meshDistributor->getRstart() * nComponents && 
+		    petscColIdx < meshDistributor->getRstart() * nComponents + nRankRows)
+		  d_nnz[localPetscRowIdx]++;
+		else
+		  o_nnz[localPetscRowIdx]++;
+	      }    
+	    }
+	  } else {
+	    // === The current row DOF is not a rank dof, i.e., it will be created ===
+	    // === on this rank, but after this it will be send to another rank    ===
+	    // === matrix. So we need to send also the corresponding nnz structure ===
+	    // === of this row to the corresponding rank.                          ===
+	    
+	    // Send all non zero entries to the member of the row DOF.
+	    int sendToRank = sendDofToRank[*cursor];
+	    
+	    for (icursor_type icursor = begin<nz>(cursor), 
+		   icend = end<nz>(cursor); icursor != icend; ++icursor) {
+	      if (value(*icursor) != 0.0) {
 		int petscColIdx = 
 		  meshDistributor->mapLocalToGlobal(col(*icursor)) * nComponents + j;
-
-		if (value(*icursor) != 0.0 || petscRowIdx == petscColIdx) {
-		  // The row DOF is a rank DOF, if also the column is a rank DOF, 
-		  // increment the d_nnz values for this row, otherwise the o_nnz value.
-		  if (petscColIdx >= meshDistributor->getRstart() * nComponents && 
-		      petscColIdx < meshDistributor->getRstart() * nComponents + nRankRows)
-		    d_nnz[localPetscRowIdx]++;
-		  else
-		    o_nnz[localPetscRowIdx]++;
-		}    
-	      }
-	    } else {
-	      typedef std::map<int, DofContainer> RankToDofContainer;
-
-	      // === The current row DOF is not a rank dof, i.e., it will be created ===
-	      // === on this rank, but after this it will be send to another rank    ===
-	      // === matrix. So we need to send also the corresponding nnz structure ===
-	      // === of this row to the corresponding rank.                          ===
-
-	      // Find out who is the member of this DOF.
-	      int sendToRank = -1;
-	      for (RankToDofContainer::iterator it = meshDistributor->getRecvDofs().begin();
-		   it != meshDistributor->getRecvDofs().end(); ++it) {
-		for (DofContainer::iterator dofIt = it->second.begin();
-		     dofIt != it->second.end(); ++dofIt) {
-		  if (**dofIt == *cursor) {
-		    sendToRank = it->first;
-		    break;
-		  }
-		}
-
-		if (sendToRank != -1)
-		  break;
-	      }
-
-	      TEST_EXIT_DBG(sendToRank != -1)("Should not happen!\n");
-
-	      // Send all non zero entries to the member of the row DOF.
-	      for (icursor_type icursor = begin<nz>(cursor), 
-		     icend = end<nz>(cursor); icursor != icend; ++icursor) {
-		if (value(*icursor) != 0.0) {
-		  int petscColIdx = 
-		    meshDistributor->mapLocalToGlobal(col(*icursor)) * nComponents + j;
-		  
-		  sendMatrixEntry[sendToRank].
-		    push_back(std::make_pair(petscRowIdx, petscColIdx));
-		}
+		
+		sendMatrixEntry[sendToRank].
+		  push_back(std::make_pair(petscRowIdx, petscColIdx));
 	      }
-
-	    } // if (isRankDof[*cursor]) ... else ...
-	  } // for each row in mat[i][j]
-	} // if mat[i][j]
+	    }
+	    
+	  } // if (isRankDof[*cursor]) ... else ...
+	} // for each row in mat[i][j]
       } 
     }
 
@@ -368,7 +400,9 @@ namespace AMDiS {
 
     StdMpi<MatrixNnzEntry> stdMpi(meshDistributor->getMpiComm(), true);
     stdMpi.send(sendMatrixEntry);
-    stdMpi.recv(meshDistributor->getSendDofs());
+    for (std::set<int>::iterator it = recvFromRank.begin(); 
+	 it != recvFromRank.end(); ++it)
+      stdMpi.recv(*it);
     stdMpi.startCommunication<int>(MPI_INT);
 
 
@@ -472,8 +506,8 @@ namespace AMDiS {
     for (int i = 0; i < nComponents; i++)
       for (int j = 0; j < nComponents; j++)
 	if ((*mat)[i][j])
-	  setDofMatrix((*mat)[i][j], nComponents, i, j);	
-	
+	  setDofMatrix((*mat)[i][j], nComponents, i, j);
+
 #if (DEBUG != 0)
     INFO(info, 8)("Fill petsc matrix 2 needed %.5f seconds\n", MPI::Wtime() - wtime);
 #endif