Fixed some problems with nnz computation with coarse grid, still not finished.

AMDiS/src/parallel/PetscSolverGlobalMatrix.cc

@@ -51,7 +51,9 @@ namespace AMDiS {
 
     MatCreateAIJ(mpiCommGlobal, nRankRows, nRankRows, 
 		 nOverallRows, nOverallRows,
-		 0, matInteriorDiagNnz, 0, matInteriorOffdiagNnz, &matIntInt);
+		 0, nnzInterior.dnnz, 
+		 0, nnzInterior.onnz, 
+		 &matIntInt);
 
     MatSetOption(matIntInt, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
 
@@ -132,28 +134,40 @@ namespace AMDiS {
 
     bool localMatrix = (subdomainLevel == 0);
     if (checkMeshChange(mat, localMatrix)) {
-      createPetscNnzStructureWithCoarseSpace(mat, 
-					     *interiorMap, 
-					     matInteriorDiagNnz,
-					     matInteriorOffdiagNnz,
-					     localMatrix);
-
-      if (coarseSpaceMap)
+       createPetscNnzStructureWithCoarseSpace(mat, 
+					      *interiorMap,
+					      *interiorMap,
+					      nnzInterior,
+					      localMatrix);
+
+      if (coarseSpaceMap) {
+ 	createPetscNnzStructureWithCoarseSpace(mat, 
+ 					       *coarseSpaceMap,
+ 					       *coarseSpaceMap,
+ 					       nnzCoarse);
 	createPetscNnzStructureWithCoarseSpace(mat, 
-					       *coarseSpaceMap, 
-					       matCoarseDiagNnz,
-					       matCoarseOffdiagNnz);	
+					       *coarseSpaceMap,
+					       *interiorMap,
+					       nnzCoarseInt);
+ 	createPetscNnzStructureWithCoarseSpace(mat, 
+ 					       *interiorMap,
+ 					       *coarseSpaceMap,
+ 					       nnzIntCoarse);
+      }
     }
 
     if (localMatrix) {
       MatCreateSeqAIJ(mpiCommLocal, nRowsRankInterior, nRowsRankInterior,
-		      0, matInteriorDiagNnz, &matIntInt);
+		      0, nnzInterior.dnnz, 
+		      &matIntInt);
+      //      MatSetOption(matIntInt, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
     } else {
       MatCreateAIJ(mpiCommLocal, 
 		   nRowsRankInterior, nRowsRankInterior,
 		   nRowsOverallInterior, nRowsOverallInterior,
-		   0, matInteriorDiagNnz, 0, matInteriorOffdiagNnz, 
+		   0, nnzInterior.dnnz, 0, nnzInterior.onnz, 
 		   &matIntInt);
+      //      MatSetOption(matIntInt, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
     }
     
     if (coarseSpaceMap) {
@@ -163,13 +177,16 @@ namespace AMDiS {
       MatCreateAIJ(mpiCommGlobal,
 		   nRowsRankCoarse, nRowsRankCoarse,
 		   nRowsOverallCoarse, nRowsOverallCoarse,
-		   0, matCoarseDiagNnz, 0, matCoarseOffdiagNnz, 
+		   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);
 
@@ -177,7 +194,7 @@ namespace AMDiS {
 		   nRowsRankInterior, nRowsRankCoarse,
 		   nGlobalOverallInterior, nRowsOverallCoarse,
 		   100, PETSC_NULL, 100, PETSC_NULL,
-		   //	   0, matInteriorDiagNnz, 0, matCoarseOffdiagNnz,
+		   //	   0, nnzIntCoarse.dnnz, 0, nnzIntCoarse.onnz,
 		   &matIntCoarse);
       MatSetOption(matIntCoarse, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
     }
@@ -259,14 +276,19 @@ namespace AMDiS {
 
 	    if (colsOther.size()) {
 	      if (subdomainLevel == 0) {
-		for (unsigned int k = 0; k < colsOther.size(); k++)
+		for (unsigned int k = 0; k < colsOther.size(); k++) {
+		  int local = interiorMap->getMatIndex(j, colsOther[k]) - interiorMap->getStartDofs();
+		  if (rowIndex == 0 && local == 9297)
+		    MSG("FOUND: %d %d %d\n", colsOther[k], j, interiorMap->getMatIndex(j, colsOther[k]));
+
 		  colsOther[k] = interiorMap->getMatIndex(j, colsOther[k]);
+		}
 	      } else {
 		for (unsigned int k = 0; k < colsOther.size(); k++)
 		  colsOther[k] = 
 		    interiorMap->getMatIndex(j, colsOther[k]) + rStartInterior;
 	      }
- 	      
+
 	      MatSetValues(matCoarseInt, 1, &rowIndex, colsOther.size(),
  			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
 	    }
@@ -318,6 +340,8 @@ namespace AMDiS {
       MatAssemblyEnd(matCoarseInt, MAT_FINAL_ASSEMBLY);
     }
 
+    AMDiS::finalize();
+    exit(0);
 
     // === Remove Dirichlet BC DOFs. ===
 
@@ -577,30 +601,18 @@ namespace AMDiS {
       static_cast<int>(meshDistributor->getLastMeshChangeIndex() != lastMeshNnz);
     mpiCommGlobal.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);
 
-    if (!matInteriorDiagNnz || recvAllValues != 0 || alwaysCreateNnzStructure) {
+    if (!nnzInterior.dnnz || recvAllValues != 0 || alwaysCreateNnzStructure) {
       vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);
       
       interiorMap->setComputeMatIndex(true, !localMatrix);
       interiorMap->update(feSpaces);
 
-      if (matInteriorDiagNnz) {
-	delete [] matInteriorDiagNnz;
-	matInteriorDiagNnz = NULL;
-      }
-      
-      if (matInteriorOffdiagNnz) {
-	delete [] matInteriorOffdiagNnz;
-	matInteriorOffdiagNnz = NULL;
-      }
+      nnzInterior.clear();
 
-      if (matCoarseDiagNnz) {
-	delete [] matCoarseDiagNnz;
-	matCoarseDiagNnz = NULL;
-      }
-      
-      if (matCoarseOffdiagNnz) {
-	delete [] matCoarseOffdiagNnz;
-	matCoarseOffdiagNnz = NULL;
+      if (coarseSpaceMap) {
+	nnzCoarse.clear();
+	nnzCoarseInt.clear();
+	nnzIntCoarse.clear();
       }
 
       updateSubdomainData();
@@ -1034,18 +1046,11 @@ namespace AMDiS {
   {
     FUNCNAME("PetscSolverGlobalMatrix::createPetscNnzStructure()");
 
-    TEST_EXIT_DBG(!matInteriorDiagNnz)("There is something wrong!\n");
-
     vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);
     int nRankRows = interiorMap->getRankDofs();
     int rankStartIndex = interiorMap->getStartDofs();
 
-    matInteriorDiagNnz = new int[nRankRows];
-    matInteriorOffdiagNnz = new int[nRankRows];
-    for (int i = 0; i < nRankRows; i++) {
-      matInteriorDiagNnz[i] = 0;
-      matInteriorOffdiagNnz[i] = 0;
-    }
+    nnzInterior.create(nRankRows, nRankRows);
 
     using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
     namespace traits = mtl::traits;
@@ -1143,9 +1148,9 @@ namespace AMDiS {
 		// otherwise the matInteriorDiagNnz value.
 		if (petscColIdx >= rankStartIndex && 
 		    petscColIdx < rankStartIndex + nRankRows)
-		  matInteriorDiagNnz[localPetscRowIdx]++;
+		  nnzInterior.dnnz[localPetscRowIdx]++;
 		else
-		  matInteriorOffdiagNnz[localPetscRowIdx]++;
+		  nnzInterior.onnz[localPetscRowIdx]++;
 	      }    
 	    }
 	  } else {
@@ -1202,9 +1207,9 @@ namespace AMDiS {
 	     r, localRowIdx, nRankRows, it->first);
 	  
 	  if (c < rankStartIndex || c >= rankStartIndex + nRankRows)
-	    matInteriorOffdiagNnz[localRowIdx]++;
+	    nnzInterior.onnz[localRowIdx]++;
 	  else
-	    matInteriorDiagNnz[localRowIdx]++;
+	    nnzInterior.dnnz[localRowIdx]++;
 	}
       }
     }
@@ -1217,33 +1222,26 @@ namespace AMDiS {
 
     if (nRankRows < 100) 
       for (int i = 0; i < nRankRows; i++)
-	matInteriorDiagNnz[i] = std::min(matInteriorDiagNnz[i], nRankRows);
+	nnzInterior.dnnz[i] = std::min(nnzInterior.dnnz[i], nRankRows);
   }
 
 
   void PetscSolverGlobalMatrix::createPetscNnzStructureWithCoarseSpace(Matrix<DOFMatrix*> *mat,
-								       ParallelDofMapping &dofMap,
-								       int*& diagNnz, 
-								       int*& offdiagNnz,
+								       ParallelDofMapping &rowDofMap,
+								       ParallelDofMapping &colDofMap,
+								       NnzStructure &nnz,
 								       bool localMatrix)
   {
     FUNCNAME("PetscSolverGlobalMatrix::createPetscNnzStructure()");
 
-    TEST_EXIT_DBG(!diagNnz)("There is something wrong!\n");
-
     vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);
-    int nRankRows = dofMap.getRankDofs();
-    int rankStartIndex = dofMap.getStartDofs();
+    int nRankRows = rowDofMap.getRankDofs();
+    int rankStartRowIndex = rowDofMap.getStartDofs();
 
-    diagNnz = new int[nRankRows];
-    for (int i = 0; i < nRankRows; i++)
-      diagNnz[i] = 0;
+    int nRankCols = colDofMap.getRankDofs();
+    int rankStartColIndex = colDofMap.getStartDofs();
 
-    if (localMatrix == false) {
-      offdiagNnz = new int[nRankRows];
-      for (int i = 0; i < nRankRows; i++)
-	offdiagNnz[i] = 0;
-    }
+    nnz.create(nRankRows, (!localMatrix ? nRankRows : -1));
 
     using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
     namespace traits = mtl::traits;
@@ -1308,21 +1306,21 @@ namespace AMDiS {
 	for (cursor_type cursor = begin<row>(bmat), 
 	       cend = end<row>(bmat); cursor != cend; ++cursor) {
 
-	  if (dofMap[feSpaces[i]].find(*cursor, rowDofIndex) == false)
+	  if (rowDofMap[feSpaces[i]].find(*cursor, rowDofIndex) == false)
 	    continue;
 
 	  // The corresponding global matrix row index of the current row DOF.
 	  int petscRowIdx = 0;
 	  if (localMatrix) {
-	    petscRowIdx = dofMap.getLocalMatIndex(i, *cursor);
+	    petscRowIdx = rowDofMap.getLocalMatIndex(i, *cursor);
 	  } else {
-	    if (dofMap.isMatIndexFromGlobal())
-	      petscRowIdx = dofMap.getMatIndex(i, rowDofIndex.global);
+	    if (rowDofMap.isMatIndexFromGlobal())
+	      petscRowIdx = rowDofMap.getMatIndex(i, rowDofIndex.global);
 	    else
-	      petscRowIdx = dofMap.getMatIndex(i, *cursor);
+	      petscRowIdx = rowDofMap.getMatIndex(i, *cursor);
 	  }
 
-	  if (localMatrix || dofMap[feSpaces[i]].isRankDof(*cursor)) {
+	  if (localMatrix || rowDofMap[feSpaces[i]].isRankDof(*cursor)) {
     	    
 	    // === The current row DOF is a rank DOF, so create the       ===
 	    // === corresponding nnz values directly on rank's nnz data.  ===
@@ -1330,50 +1328,47 @@ namespace AMDiS {
 	    // This is the local row index of the local PETSc matrix.
 	    int localPetscRowIdx = petscRowIdx;
 
-	    if (localMatrix == false) {
-	      localPetscRowIdx -= rankStartIndex;
+	    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,
-		 dofMap[feSpaces[i]][*cursor].global,
-		 petscRowIdx, 
-		 localPetscRowIdx, 
-		 rankStartIndex,
-		 i,
-		 nComponents, 
-		 nRankRows);
-	    }
+	    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[feSpaces[i]][*cursor].global,
+	       petscRowIdx, 
+	       localPetscRowIdx, 
+	       rankStartRowIndex,
+	       i,
+	       nComponents, 
+	       nRankRows);
 	    
 
 	    if (localMatrix) {
 	      for (icursor_type icursor = begin<nz>(cursor), 
 		     icend = end<nz>(cursor); icursor != icend; ++icursor)
-		if (dofMap[feSpaces[j]].isSet(col(*icursor)))
-		  diagNnz[localPetscRowIdx]++;	      
+		if (colDofMap[feSpaces[j]].isSet(col(*icursor)))
+		  nnz.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 (dofMap[feSpaces[j]].find(col(*icursor), colDofIndex) == false)
+		if (colDofMap[feSpaces[j]].find(col(*icursor), colDofIndex) == false)
 		  continue;
 		
-		int petscColIdx = (dofMap.isMatIndexFromGlobal() ? 
-				   dofMap.getMatIndex(j, colDofIndex.global) :
-				   dofMap.getMatIndex(j, col(*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 diagNnz values for this row, 
-		  // otherwise the offdiagNnz value.
-		  if (petscColIdx >= rankStartIndex && 
-		      petscColIdx < rankStartIndex + nRankRows)
-		    diagNnz[localPetscRowIdx]++;
-		  else
-		    offdiagNnz[localPetscRowIdx]++;
-		}    
+		int petscColIdx = (colDofMap.isMatIndexFromGlobal() ? 
+				   colDofMap.getMatIndex(j, colDofIndex.global) :
+				   colDofMap.getMatIndex(j, 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)
+		  nnz.dnnz[localPetscRowIdx]++;
+		else
+		  nnz.onnz[localPetscRowIdx]++;
 	      }
 	    }
 	  } else {
@@ -1389,17 +1384,15 @@ namespace AMDiS {
 
 	    for (icursor_type icursor = begin<nz>(cursor), 
 		   icend = end<nz>(cursor); icursor != icend; ++icursor) {
-	      if (dofMap[feSpaces[j]].find(col(*icursor), colDofIndex) == false)
+	      if (colDofMap[feSpaces[j]].find(col(*icursor), colDofIndex) == false)
 		continue;
 
-	      if (value(*icursor) != 0.0) {
-		int petscColIdx = (dofMap.isMatIndexFromGlobal() ? 
-				   dofMap.getMatIndex(j, colDofIndex.global) :
-				   dofMap.getMatIndex(j, col(*icursor)));
-		
-		sendMatrixEntry[sendToRank].
-		  push_back(make_pair(petscRowIdx, petscColIdx));
-	      }
+	      int petscColIdx = (colDofMap.isMatIndexFromGlobal() ? 
+				 colDofMap.getMatIndex(j, colDofIndex.global) :
+				 colDofMap.getMatIndex(j, col(*icursor)));
+	      
+	      sendMatrixEntry[sendToRank].
+		push_back(make_pair(petscRowIdx, petscColIdx));
 	    }
 	    
 	  } // if (isRankDof[*cursor]) ... else ...
@@ -1407,6 +1400,7 @@ namespace AMDiS {
       } 
     }
 
+
     if (localMatrix == false) {
       // === Send and recv the nnz row structure to/from other ranks. ===
       
@@ -1428,43 +1422,54 @@ namespace AMDiS {
 	    int r = it->second[i].first;
 	    int c = it->second[i].second;
 	    
-	    int localRowIdx = r - rankStartIndex;
+	    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 < rankStartIndex || c >= rankStartIndex + nRankRows)
-	      offdiagNnz[localRowIdx]++;
+	    if (c < rankStartColIndex || c >= rankStartColIndex + nRankCols)
+	      nnz.onnz[localRowIdx]++;
 	    else
-	      diagNnz[localRowIdx]++;
+	      nnz.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++)
-	  diagNnz[i] = std::min(diagNnz[i], nRankRows);
     }
 
+    // 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], nRankCols);
 
-#if (DEBUG != 0)    
+#if (DEBUG != 0)
     int nMax = 0; 
     int nSum = 0;
     for (int i = 0; i < nRankRows; i++) {
-      nMax = std::max(nMax, diagNnz[i]);
-      nSum += diagNnz[i];
+      nMax = std::max(nMax, nnz.dnnz[i]);
+      nSum += nnz.dnnz[i];
     }
 
-    MSG("NNZ in matrix: max = %d, avrg = %.0f\n",
-	nMax, static_cast<double>(nSum) / nRankRows);
-#endif
+    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
 }

AMDiS/src/parallel/PetscSolverGlobalMatrix.h

@@ -31,16 +31,56 @@ namespace AMDiS {
 
   using namespace std;
 
+  struct NnzStructure {
+    NnzStructure()
+      : dnnz(NULL),