From 2bbf10f76aec8b98ce655056ee4a22d7de296059 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Tue, 26 Jun 2012 20:11:48 +0000
Subject: [PATCH] Fixed some problems with nnz computation with coarse grid,
still not finished.
---
AMDiS/src/parallel/PetscSolverGlobalMatrix.cc | 261 +++++++++---------
AMDiS/src/parallel/PetscSolverGlobalMatrix.h | 62 ++++-
2 files changed, 184 insertions(+), 139 deletions(-)
diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
index 06552c19..68033d96 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
+++ b/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
}
diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.h b/AMDiS/src/parallel/PetscSolverGlobalMatrix.h
index 63b2fc47..728b1a2e 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.h
+++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.h
@@ -31,16 +31,56 @@ namespace AMDiS {
using namespace std;
+ struct NnzStructure {
+ NnzStructure()
+ : dnnz(NULL),
+ onnz(NULL)
+ {}
+
+ void clear()
+ {
+ if (dnnz) {
+ delete [] dnnz;
+ dnnz = NULL;
+ }
+
+ if (onnz) {
+ delete [] onnz;
+ onnz = NULL;
+ }
+ }
+
+ void create(int nRows0, int nRows1 = -1)
+ {
+ 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;
+ }
+ }
+
+ int *dnnz;
+
+ int *onnz;
+ };
+
+
+
class PetscSolverGlobalMatrix : public PetscSolver
{
public:
PetscSolverGlobalMatrix()
: PetscSolver(),
petscSolVec(PETSC_NULL),
- matInteriorDiagNnz(NULL),
- matInteriorOffdiagNnz(NULL),
- matCoarseDiagNnz(NULL),
- matCoarseOffdiagNnz(NULL),
lastMeshNnz(0),
zeroStartVector(false),
alwaysCreateNnzStructure(false)
@@ -81,9 +121,9 @@ namespace AMDiS {
void createPetscNnzStructure(Matrix<DOFMatrix*> *mat);
void createPetscNnzStructureWithCoarseSpace(Matrix<DOFMatrix*> *mat,
- ParallelDofMapping &dofMap,
- int*& diagNnz,
- int*& offdiagNnz,
+ ParallelDofMapping &rowDofMap,
+ ParallelDofMapping &colDofMap,
+ NnzStructure &nnz,
bool localMatrix = false);
/// Takes a DOF matrix and sends the values to the global PETSc matrix.
@@ -108,11 +148,11 @@ namespace AMDiS {
protected:
- /// Arrays definig the non zero pattern of Petsc's interior matrix.
- int *matInteriorDiagNnz, *matInteriorOffdiagNnz;
+ NnzStructure nnzInterior;
+
+ NnzStructure nnzCoarse;
- /// Arrays definig the non zero pattern of Petsc's coarse space matrix.
- int *matCoarseDiagNnz, *matCoarseOffdiagNnz;
+ NnzStructure nnzIntCoarse, nnzCoarseInt;
Vec petscSolVec;
--
GitLab