Merge.

AMDiS/src/parallel/MatrixNnzStructure.cc

@@ -10,6 +10,7 @@
 // See also license.opensource.txt in the distribution.
 
 
+#include <algorithm>
 #include "Global.h"
 #include "DOFMatrix.h"
 #include "parallel/MatrixNnzStructure.h"
@@ -286,8 +287,8 @@ namespace AMDiS {
       nSum = 0;
 
       for (int i = 0; i < nRankRows; i++) {
-	nMax = std::max(nMax, nnz.onnz[i]);
-	nSum += nnz.onnz[i];
+	nMax = std::max(nMax, onnz[i]);
+	nSum += onnz[i];
       }
 
       MSG("NNZ in offdiag block: max = %d, avrg = %.0f\n",

AMDiS/src/parallel/PetscSolverGlobalMatrix.cc

@@ -1016,189 +1016,4 @@ namespace AMDiS {
     }
   }
 
-#if 0
-  void PetscSolverGlobalMatrix::createPetscNnzStructure(Matrix<DOFMatrix*> *mat)
-  {
-    FUNCNAME("PetscSolverGlobalMatrix::createPetscNnzStructure()");
-
-    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);
-    int nRankRows = interiorMap->getRankDofs();
-    int rankStartIndex = interiorMap->getStartDofs();
-
-    nnzInterior.create(nRankRows, nRankRows);
-
-    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;
-	
-	for (cursor_type cursor = begin<row>(bmat), 
-	       cend = end<row>(bmat); cursor != cend; ++cursor) {
-	  int globalRowDof = (*interiorMap)[feSpaces[i]][*cursor].global;
-
-	  // The corresponding global matrix row index of the current row DOF.
-	  int petscRowIdx = interiorMap->getMatIndex(i, globalRowDof);
-	  if ((*interiorMap)[feSpaces[i]].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 - rankStartIndex;
-	    
-	    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,
-	       (*interiorMap)[feSpaces[i]][*cursor].global,
-	       petscRowIdx, 
-	       localPetscRowIdx, 
-	       rankStartIndex,
-	       i,
-	       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 globalColDof = (*interiorMap)[feSpaces[j]][col(*icursor)].global;
-	      int petscColIdx = interiorMap->getMatIndex(j, globalColDof);
-	      
-	      if (value(*icursor) != 0.0 || petscRowIdx == petscColIdx) {
-		// The row DOF is a rank DOF, if also the column is a rank DOF, 
-		// increment the matInteriorDiagNnz values for this row, 
-		// otherwise the matInteriorDiagNnz value.
-		if (petscColIdx >= rankStartIndex && 
-		    petscColIdx < rankStartIndex + nRankRows)
-		  nnzInterior.dnnz[localPetscRowIdx]++;
-		else
-		  nnzInterior.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, i)];
-	    
-	    for (icursor_type icursor = begin<nz>(cursor), 
-		   icend = end<nz>(cursor); icursor != icend; ++icursor) {
-	      if (value(*icursor) != 0.0) {
-		int globalColDof = 
-		  (*interiorMap)[feSpaces[j]][col(*icursor)].global;
-		int petscColIdx = interiorMap->getMatIndex(j, globalColDof);
-		
-		sendMatrixEntry[sendToRank].
-		  push_back(make_pair(petscRowIdx, petscColIdx));
-	      }
-	    }
-	    
-	  } // if (isRankDof[*cursor]) ... else ...
-	} // for each row in mat[i][j]
-      } 
-    }
-
-    // === Send and recv the nnz row structure to/from other ranks. ===
-
-    StdMpi<MatrixNnzEntry> stdMpi(mpiCommGlobal, 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 - rankStartIndex;
-
-	  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 < rankStartIndex || c >= rankStartIndex + nRankRows)
-	    nnzInterior.onnz[localRowIdx]++;
-	  else
-	    nnzInterior.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++)
-	nnzInterior.dnnz[i] = std::min(nnzInterior.dnnz[i], nRankRows);
-  }
-#endif
-
 }

AMDiS/src/parallel/PetscSolverGlobalMatrix.h

@@ -75,9 +75,6 @@ namespace AMDiS {
 
     virtual void exitPreconditioner(PC pc);
 
-    /// Creates a new non zero pattern structure for the PETSc matrix. 
-    //    void createPetscNnzStructure(Matrix<DOFMatrix*> *mat);
-
     /// Takes a DOF matrix and sends the values to the global PETSc matrix.
     void setDofMatrix(DOFMatrix* mat, int nRowMat = 0, int nColMat = 0);