Added new class to create nnz structure of a DOF matrix.

AMDiS/CMakeLists.txt

@@ -225,6 +225,7 @@ if(ENABLE_PARALLEL_DOMAIN)
 		${SOURCE_DIR}/parallel/CheckerPartitioner.cc
 		${SOURCE_DIR}/parallel/ElementObjectDatabase.cc
 		${SOURCE_DIR}/parallel/InteriorBoundary.cc
+		${SOURCE_DIR}/parallel/MatrixNnzStructure.cc
 		${SOURCE_DIR}/parallel/MeshDistributor.cc 
 		${SOURCE_DIR}/parallel/MeshLevelData.cc 
 		${SOURCE_DIR}/parallel/MeshManipulation.cc

AMDiS/src/parallel/MatrixNnzStructure.cc

+//
+// Software License for AMDiS
+//
+// Copyright (c) 2010 Dresden University of Technology 
+// All rights reserved.
+// Authors: Simon Vey, Thomas Witkowski et al.
+//
+// This file is part of AMDiS
+//
+// See also license.opensource.txt in the distribution.
+
+
+#include "Global.h"
+#include "DOFMatrix.h"
+#include "parallel/MatrixNnzStructure.h"
+#include "parallel/ParallelDofMapping.h"
+
+namespace AMDiS {
+
+  void MatrixNnzStructure::clear()
+  {
+    if (dnnz) {
+      delete [] dnnz;
+      dnnz = NULL;
+    }
+    
+    if (onnz) {
+      delete [] onnz;
+      onnz = NULL;
+    }
+  }
+
+
+  void MatrixNnzStructure::create(Matrix<DOFMatrix*> *mat,
+				  MPI::Intracomm mpiComm,
+				  ParallelDofMapping &rowDofMap,
+				  ParallelDofMapping &colDofMap,
+				  bool localMatrix)
+  {
+    FUNCNAME("MatrixNnzStructure::create()");
+
+    int nRankRows = rowDofMap.getRankDofs();
+    int rankStartRowIndex = rowDofMap.getStartDofs();
+
+    int nRankCols = colDofMap.getRankDofs();
+    int rankStartColIndex = colDofMap.getStartDofs();
+
+    create(nRankRows, (!localMatrix ? nRankRows : -1));
+
+    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
+    namespace traits = mtl::traits;
+    typedef DOFMatrix::base_matrix_type Matrix;
+    typedef traits::range_generator<row, Matrix>::type cursor_type;
+    typedef traits::range_generator<nz, cursor_type>::type icursor_type;
+
+    typedef vector<pair<int, int> > MatrixNnzEntry;
+    typedef 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. 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.
+    map<int, MatrixNnzEntry> sendMatrixEntry;
+
+    // Maps to each DOF that must be send to another rank the rank number of the
+    // receiving rank.
+    map<pair<DegreeOfFreedom, int>, int> sendDofToRank;
+
+    int nComponents = mat->getNumRows();
+
+    // First, create for all ranks, to which we send data to, MatrixNnzEntry 
+    // object with 0 entries.
+    for (int i = 0; i < nComponents; i++) {
+      const FiniteElemSpace* feSpace = NULL;
+      for (int j = 0; j < nComponents; j++)
+	if ((*mat)[i][j])
+	  feSpace = (*mat)[i][j]->getRowFeSpace();
+
+      for (DofComm::Iterator it(rowDofMap.getDofComm().getRecvDofs(), feSpace);
+	   !it.end(); it.nextRank()) {
+	sendMatrixEntry[it.getRank()].resize(0);
+	
+	for (; !it.endDofIter(); it.nextDof())
+	  sendDofToRank[make_pair(it.getDofIndex(), i)] = it.getRank();
+      }
+    }
+
+    // Create list of ranks from which we receive data from.
+    std::set<int> recvFromRank;
+    for (int i = 0; i < nComponents; i++)  {
+      const FiniteElemSpace* feSpace = NULL;
+      for (int j = 0; j < nComponents; j++)
+	if ((*mat)[i][j])
+	  feSpace = (*mat)[i][j]->getRowFeSpace();
+
+      for (DofComm::Iterator it(rowDofMap.getDofComm().getSendDofs(), feSpace);
+	   !it.end(); it.nextRank())
+	recvFromRank.insert(it.getRank());
+    }
+
+
+    // === Traverse matrices to create nnz data. ===
+
+    for (int rowComp = 0; rowComp < nComponents; rowComp++) {
+      for (int colComp = 0; colComp < nComponents; colComp++) {
+	DOFMatrix *dofMat = (*mat)[rowComp][colComp];
+
+ 	if (!dofMat)
+	  continue;
+
+	const FiniteElemSpace *rowFeSpace = dofMat->getRowFeSpace();
+	const FiniteElemSpace *colFeSpace = dofMat->getColFeSpace();
+
+
+	// === Prepare MTL4 iterators for the current DOF matrix. ===
+
+	Matrix baseMat = dofMat->getBaseMatrix();
+	traits::col<Matrix>::type col(baseMat);
+	traits::const_value<Matrix>::type value(baseMat);  
+	cursor_type cursor = begin<row>(baseMat);
+	cursor_type cend = end<row>(baseMat);
+
+
+	// === Iterate on all DOFs of the row mapping. ===
+
+	DofMap::iterator rowIt = rowDofMap[rowFeSpace].begin();
+	DofMap::iterator rowEndIt = rowDofMap[rowFeSpace].end();
+	for (; rowIt != rowEndIt; ++rowIt) {
+
+	  // Go to the corresponding matrix row (note, both the mapping and the
+	  // matrix rows are stored in increasing order).
+	  while (*cursor != rowIt->first)
+	    ++cursor;
+
+	  // The corresponding global matrix row index of the current row DOF.
+	  int petscRowIdx = 0;
+	  if (localMatrix) {
+	    petscRowIdx = rowDofMap.getLocalMatIndex(rowComp, *cursor);
+	  } else {
+	    if (rowDofMap.isMatIndexFromGlobal())
+	      petscRowIdx = rowDofMap.getMatIndex(rowComp, rowIt->second.global);
+	    else
+	      petscRowIdx = rowDofMap.getMatIndex(rowComp, *cursor);
+	  }
+
+
+	  if (localMatrix || rowDofMap[rowFeSpace].isRankDof(*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;
+
+	    if (localMatrix == false)
+	      localPetscRowIdx -= rankStartRowIndex;
+	    
+	    TEST_EXIT_DBG(localPetscRowIdx >= 0 && localPetscRowIdx < nRankRows)
+	      ("Should not happen! \n Debug info: DOF = %d   globalRowIndx = %d   petscRowIdx = %d   localPetscRowIdx = %d   rStart = %d   compontens = %d from %d   nRankRows = %d\n",
+	       *cursor,
+	       rowDofMap[rowFeSpace][*cursor].global,
+	       petscRowIdx, 
+	       localPetscRowIdx, 
+	       rankStartRowIndex,
+	       rowComp,
+	       nComponents, 
+	       nRankRows);
+	    
+
+	    if (localMatrix) {
+	      for (icursor_type icursor = begin<nz>(cursor), 
+		     icend = end<nz>(cursor); icursor != icend; ++icursor)
+		if (colDofMap[colFeSpace].isSet(col(*icursor)))
+		  dnnz[localPetscRowIdx]++;	      
+	    } else {
+	      MultiIndex colDofIndex;
+
+	      // Traverse all non zero entries in this row.
+	      for (icursor_type icursor = begin<nz>(cursor), 
+		     icend = end<nz>(cursor); icursor != icend; ++icursor) {
+		if (colDofMap[colFeSpace].find(col(*icursor), colDofIndex) == false)
+		  continue;
+		
+		int petscColIdx = (colDofMap.isMatIndexFromGlobal() ? 
+				   colDofMap.getMatIndex(colComp, colDofIndex.global) :
+				   colDofMap.getMatIndex(colComp, col(*icursor)));
+	
+		// The row DOF is a rank DOF, if also the column is a rank DOF, 
+		// increment the diagNnz values for this row, 
+		// otherwise the offdiagNnz value.
+		if (petscColIdx >= rankStartColIndex && 
+		    petscColIdx < rankStartColIndex + nRankCols)
+		  dnnz[localPetscRowIdx]++;
+		else
+		  onnz[localPetscRowIdx]++;
+	      }
+	    }
+	  } else {
+	    // === The current row DOF is not a rank DOF, i.e., its values   ===
+	    // === are also created on this rank, but afterthere they will   ===
+	    // === be send to another rank. So we need to send also the      ===
+	    // === corresponding nnz structure of this row to the corres-    ===
+	    // === ponding rank.                                             ===
+	    
+	    // Send all non zero entries to the member of the row DOF.
+	    int sendToRank = sendDofToRank[make_pair(*cursor, rowComp)];
+	    MultiIndex colDofIndex;
+
+	    for (icursor_type icursor = begin<nz>(cursor), 
+		   icend = end<nz>(cursor); icursor != icend; ++icursor) {
+	      if (colDofMap[colFeSpace].find(col(*icursor), colDofIndex) == false)
+		continue;
+
+	      int petscColIdx = (colDofMap.isMatIndexFromGlobal() ? 
+				 colDofMap.getMatIndex(colComp, colDofIndex.global) :
+				 colDofMap.getMatIndex(colComp, col(*icursor)));
+	      
+	      sendMatrixEntry[sendToRank].
+		push_back(make_pair(petscRowIdx, petscColIdx));
+	    }
+	    
+	  } // if (isRankDof[*cursor]) ... else ...
+	} // for each row in mat[i][j]
+      } 
+    }
+
+
+    if (localMatrix == false) {
+      // === Send and recv the nnz row structure to/from other ranks. ===
+      
+      StdMpi<MatrixNnzEntry> stdMpi(mpiComm, true);
+      stdMpi.send(sendMatrixEntry);
+      for (std::set<int>::iterator it = recvFromRank.begin(); 
+	   it != recvFromRank.end(); ++it)
+	stdMpi.recv(*it);
+      stdMpi.startCommunication();
+      
+      
+      // === Evaluate the nnz structure this rank got from other ranks and add ===
+      // === it to the PETSc nnz data structure.                               ===
+      
+      for (map<int, MatrixNnzEntry>::iterator it = stdMpi.getRecvData().begin();
+	   it != stdMpi.getRecvData().end(); ++it) {
+	if (it->second.size() > 0) {
+	  for (unsigned int i = 0; i < it->second.size(); i++) {
+	    int r = it->second[i].first;
+	    int c = it->second[i].second;
+	    
+	    int localRowIdx = r - rankStartRowIndex;
+	    
+	    TEST_EXIT_DBG(localRowIdx >= 0 && localRowIdx < nRankRows)
+	      ("Got row index %d/%d (nRankRows = %d) from rank %d. Should not happen!\n",
+	       r, localRowIdx, nRankRows, it->first);
+	    
+	    if (c < rankStartColIndex || c >= rankStartColIndex + nRankCols)
+	      onnz[localRowIdx]++;
+	    else
+	      dnnz[localRowIdx]++;
+	  }
+	}
+      }
+    }
+
+    // The above algorithm for calculating the number of nnz per row over-
+    // approximates the value, i.e., the number is always equal or larger to 
+    // the real number of nnz values in the global parallel matrix. For small
+    // matrices, the problem may arise, that the result is larger than the
+    // number of elements in a row. This is fixed in the following.
+
+    if (nRankRows < 100) 
+      for (int i = 0; i < nRankRows; i++)
+ 	dnnz[i] = std::min(dnnz[i], nRankCols);
+
+#if (DEBUG != 0)
+    int nMax = 0; 
+    int nSum = 0;
+    for (int i = 0; i < nRankRows; i++) {
+      nMax = std::max(nMax, dnnz[i]);
+      nSum += dnnz[i];
+    }
+
+    MSG("NNZ in diag block: max = %d, avrg = %.0f\n",
+	nMax, (nSum > 0 ? (static_cast<double>(nSum) / nRankRows) : 0));
+
+    if (!localMatrix) {
+      nMax = 0;
+      nSum = 0;
+
+      for (int i = 0; i < nRankRows; i++) {
+	nMax = std::max(nMax, nnz.onnz[i]);
+	nSum += nnz.onnz[i];
+      }
+
+      MSG("NNZ in offdiag block: max = %d, avrg = %.0f\n",
+	  nMax, (nSum > 0 ? (static_cast<double>(nSum) / nRankRows) : 0));     
+    }
+#endif
+  }
+
+
+  void MatrixNnzStructure::create(int nRows0, int nRows1)
+  {
+    if (nRows0 == 0)
+      return;
+    
+    TEST_EXIT_DBG(nRows0 > 0)("Should not happen!\n");
+    
+    dnnz = new int[nRows0];
+    for (int i = 0; i < nRows0; i++)
+      dnnz[i] = 0;
+    
+    if (nRows1 > 0) {
+      onnz = new int[nRows1];
+      for (int i = 0; i < nRows1; i++)
+	onnz[i] = 0;	
+    }
+  }
+
+}

AMDiS/src/parallel/MatrixNnzStructure.h

+// ============================================================================
+// ==                                                                        ==
+// == AMDiS - Adaptive multidimensional simulations                          ==
+// ==                                                                        ==
+// ==  http://www.amdis-fem.org                                              ==
+// ==                                                                        ==
+// ============================================================================
+//
+// Software License for AMDiS
+//
+// Copyright (c) 2010 Dresden University of Technology 
+// All rights reserved.
+// Authors: Simon Vey, Thomas Witkowski et al.
+//
+// This file is part of AMDiS
+//
+// See also license.opensource.txt in the distribution.
+
+
+
+/** \file MatrixNnzStructure.h */
+
+#ifndef AMDIS_MATRIX_NNZ_STRUCTURE_H
+#define AMDIS_MATRIX_NNZ_STRUCTURE_H
+
+#include "AMDiS_fwd.h"
+#include "parallel/ParallelDofMapping.h"
+
+namespace AMDiS {
+
+  using namespace std;
+
+  class MatrixNnzStructure {
+  public:
+    MatrixNnzStructure()
+      : dnnz(NULL),
+	onnz(NULL)
+    {}
+	
+    void clear();
+
+    void create(Matrix<DOFMatrix*> *mat,
+		MPI::Intracomm mpiComm,
+		ParallelDofMapping &rowDofMap,
+		ParallelDofMapping &colDofMap,
+		bool localMatrix = false);
+
+    void create(Matrix<DOFMatrix*> *mat,
+		MPI::Intracomm mpiComm,
+		ParallelDofMapping &dofMap,
+		bool localMatrix = false)
+    {
+      create(mat, mpiComm, dofMap, dofMap, localMatrix);
+    }
+
+  protected:
+    void create(int nRows0, int nRows1 = -1);
+
+  public:
+    int *dnnz;
+
+    int *onnz;
+  };
+
+}
+
+#endif

AMDiS/src/parallel/ParallelDofMapping.h

@@ -204,6 +204,16 @@ namespace AMDiS {
       return dofMap;
     }
 
+    DofMap::iterator begin()
+    {
+      return dofMap.begin();
+    }
+
+    DofMap::iterator end()
+    {
+      return dofMap.end();
+    }
+
     /// Recomputes the mapping.
     void update();
 
@@ -226,7 +236,7 @@ namespace AMDiS {
       needGlobalMapping = b;
     }
 
-    /// Sets the DOF communicators.
+    /// Sets the DOF communicator.
     void setDofComm(DofComm &dc)
     {
       dofComm = &dc;
@@ -333,6 +343,14 @@ namespace AMDiS {
     /// about DOFs that are on interior boundaries.
     void setDofComm(DofComm &dofComm);
 
+    /// Returns the DOF communicator.
+    DofComm& getDofComm()
+    {
+      TEST_EXIT_DBG(dofComm);
+
+      return *dofComm;
+    }
+
     void setComputeMatIndex(bool b, bool global = false)
     {
       needMatIndex = b;

AMDiS/src/parallel/PetscSolverGlobalMatrix.cc

@@ -36,8 +36,7 @@ namespace AMDiS {
     // === If required, recompute non zero structure of the matrix. ===
 
     if (checkMeshChange(mat))
-      createPetscNnzStructure(mat);
-
+      nnzInterior.create(mat, mpiCommGlobal, *interiorMap);
 
     // === Create PETSc vector (solution and a temporary vector). ===
 
@@ -134,26 +133,12 @@ namespace AMDiS {
 
     bool localMatrix = (subdomainLevel == 0);
     if (checkMeshChange(mat, localMatrix)) {
-       createPetscNnzStructureWithCoarseSpace(mat, 
-					      *interiorMap,
-					      *interiorMap,
-					      nnzInterior,
-					      localMatrix);
+      nnzInterior.create(mat, mpiCommGlobal, *interiorMap, localMatrix);
 
       if (coarseSpaceMap) {
- 	createPetscNnzStructureWithCoarseSpace(mat, 
- 					       *coarseSpaceMap,
- 					       *coarseSpaceMap,
- 					       nnzCoarse);
-
-	createPetscNnzStructureWithCoarseSpace(mat, 
-					       *coarseSpaceMap,
-					       *interiorMap,
-					       nnzCoarseInt);
- 	createPetscNnzStructureWithCoarseSpace(mat, 
- 					       *interiorMap,
- 					       *coarseSpaceMap,
- 					       nnzIntCoarse);
+	nnzCoarse.create(mat, mpiCommGlobal, *coarseSpaceMap);
+	nnzCoarseInt.create(mat, mpiCommGlobal, *coarseSpaceMap, *interiorMap);
+	nnzIntCoarse.create(mat, mpiCommGlobal, *interiorMap, *coarseSpaceMap);
       }
     }
 
@@ -161,14 +146,12 @@ namespace AMDiS {
       MatCreateSeqAIJ(mpiCommLocal, nRowsRankInterior, nRowsRankInterior,
 		      0, nnzInterior.dnnz, 
 		      &matIntInt);
-      //      MatSetOption(matIntInt, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
     } else {
       MatCreateAIJ(mpiCommLocal, 
 		   nRowsRankInterior, nRowsRankInterior,
 		   nRowsOverallInterior, nRowsOverallInterior,
 		   0, nnzInterior.dnnz, 0, nnzInterior.onnz, 
 		   &matIntInt);
-      //      MatSetOption(matIntInt, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
     }
     
     if (coarseSpaceMap) {
@@ -179,25 +162,19 @@ namespace AMDiS {
 		   nRowsRankCoarse, nRowsRankCoarse,
 		   nRowsOverallCoarse, nRowsOverallCoarse,
 		   0, nnzCoarse.dnnz, 0, nnzCoarse.onnz, 
-		   //	   100, PETSC_NULL, 100, PETSC_NULL,
 		   &matCoarseCoarse);
-      //      MatSetOption(matCoarseCoarse, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
 
       MatCreateAIJ(mpiCommGlobal,
 		   nRowsRankCoarse, nRowsRankInterior,
 		   nRowsOverallCoarse, nGlobalOverallInterior,
-		   //		   100, PETSC_NULL, 100, PETSC_NULL,
 		   0, nnzCoarseInt.dnnz, 0, nnzCoarseInt.onnz,
 		   &matCoarseInt);
-      //      MatSetOption(matCoarseInt, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
 
       MatCreateAIJ(mpiCommGlobal,
 		   nRowsRankInterior, nRowsRankCoarse,
 		   nGlobalOverallInterior, nRowsOverallCoarse,
-		   //		   100, PETSC_NULL, 100, PETSC_NULL,
 		   0, nnzIntCoarse.dnnz, 0, nnzIntCoarse.onnz,
 		   &matIntCoarse);
-      //      MatSetOption(matIntCoarse, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
     }
 
     // === Prepare traverse of sequentially created matrices. ===
@@ -1039,7 +1016,7 @@ namespace AMDiS {
     }
   }
 
-
+#if 0
   void PetscSolverGlobalMatrix::createPetscNnzStructure(Matrix<DOFMatrix*> *mat)
   {
     FUNCNAME("PetscSolverGlobalMatrix::createPetscNnzStructure()");
@@ -1222,253 +1199,6 @@ namespace AMDiS {
       for (int i = 0; i < nRankRows; i++)
 	nnzInterior.dnnz[i] = std::min(nnzInterior.dnnz[i], nRankRows);
   }
-
-
-  void PetscSolverGlobalMatrix::createPetscNnzStructureWithCoarseSpace(Matrix<DOFMatrix*> *mat,
-								       ParallelDofMapping &rowDofMap,
-								       ParallelDofMapping &colDofMap,
-								       NnzStructure &nnz,
-								       bool localMatrix)
-  {
-    FUNCNAME("PetscSolverGlobalMatrix::createPetscNnzStructure()");
-
-    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);
-    int nRankRows = rowDofMap.getRankDofs();
-    int rankStartRowIndex = rowDofMap.getStartDofs();
-
-    int nRankCols = colDofMap.getRankDofs();
-    int rankStartColIndex = colDofMap.getStartDofs();
-
-    nnz.create(nRankRows, (!localMatrix ? nRankRows : -1));
-
-    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
-    namespace traits = mtl::traits;
-    typedef DOFMatrix::base_matrix_type Matrix;
-    typedef vector<pair<int, int> > MatrixNnzEntry;
-    typedef 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. 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.
-    map<int, MatrixNnzEntry> sendMatrixEntry;
-
-    // Maps to each DOF that must be send to another rank the rank number of the
-    // receiving rank.
-    map<pair<DegreeOfFreedom, int>, int> sendDofToRank;
-
-
-    // First, create for all ranks, to which we send data to, MatrixNnzEntry 
-    // object with 0 entries.
-    for (unsigned int i = 0; i < feSpaces.size(); i++) {
-      for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), feSpaces[i]);
-	   !it.end(); it.nextRank()) {
-	sendMatrixEntry[it.getRank()].resize(0);
-	
-	for (; !it.endDofIter(); it.nextDof())
-	  sendDofToRank[make_pair(it.getDofIndex(), i)] = it.getRank();
-      }
-    }
-
-    // Create list of ranks from which we receive data from.
-    std::set<int> recvFromRank;
-    for (unsigned int i = 0; i < feSpaces.size(); i++) 
-      for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), feSpaces[i]);
-	   !it.end(); it.nextRank())
-	recvFromRank.insert(it.getRank());
-
-
-    // === Traverse matrices to create nnz data. ===
-
-    int nComponents = mat->getNumRows();
-    for (int i = 0; i < nComponents; i++) {
-      for (int j = 0; j < nComponents; j++) {
- 	if (!(*mat)[i][j])
-	  continue;
-
-	TEST_EXIT_DBG((*mat)[i][j]->getRowFeSpace() == feSpaces[i])
-	  ("Should not happen!\n");
-	TEST_EXIT_DBG((*mat)[i][j]->getColFeSpace() == feSpaces[j])
-	  ("Should not happen!\n");
-
-	Matrix bmat = (*mat)[i][j]->getBaseMatrix();
-
-	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;