diff --git a/AMDiS/CMakeLists.txt b/AMDiS/CMakeLists.txt
index 4e91547610b48b4f0f262ff984a2512496b7eca3..3b9b8190a3b362d1d44b7f1da6d040ad3699c4d3 100644
--- a/AMDiS/CMakeLists.txt
+++ b/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
diff --git a/AMDiS/src/parallel/MatrixNnzStructure.cc b/AMDiS/src/parallel/MatrixNnzStructure.cc
new file mode 100644
index 0000000000000000000000000000000000000000..82cc355a3c3036e173c54f403ec4583d97441ab5
--- /dev/null
+++ b/AMDiS/src/parallel/MatrixNnzStructure.cc
@@ -0,0 +1,318 @@
+//
+// 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;
+ }
+ }
+
+}
diff --git a/AMDiS/src/parallel/MatrixNnzStructure.h b/AMDiS/src/parallel/MatrixNnzStructure.h
new file mode 100644
index 0000000000000000000000000000000000000000..ffce7113a1b199fc75d82047a09d21d662653b2e
--- /dev/null
+++ b/AMDiS/src/parallel/MatrixNnzStructure.h
@@ -0,0 +1,67 @@
+// ============================================================================
+// == ==
+// == 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
diff --git a/AMDiS/src/parallel/ParallelDofMapping.h b/AMDiS/src/parallel/ParallelDofMapping.h
index 32eea63035ecc665ffd3f83b0b883756f5828895..7da23a1a31c5e2fedee5248d498417d9af96a2eb 100644
--- a/AMDiS/src/parallel/ParallelDofMapping.h
+++ b/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;
diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
index 3b547ba236bae42ddc89811bf0fe20fa4c2dbe36..418fb714a173cc9c3afcb55d1cc46f5fd46ea30e 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
+++ 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;
-
- MultiIndex rowDofIndex;
-
- for (cursor_type cursor = begin<row>(bmat),
- cend = end<row>(bmat); cursor != cend; ++cursor) {
-
- 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 = rowDofMap.getLocalMatIndex(i, *cursor);
- } else {
- if (rowDofMap.isMatIndexFromGlobal())
- petscRowIdx = rowDofMap.getMatIndex(i, rowDofIndex.global);
- else
- petscRowIdx = rowDofMap.getMatIndex(i, *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. ===
-
- // 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[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 (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 (colDofMap[feSpaces[j]].find(col(*icursor), colDofIndex) == false)
- continue;
-
- 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 {
- // === 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)];
- MultiIndex colDofIndex;
-
- for (icursor_type icursor = begin<nz>(cursor),
- icend = end<nz>(cursor); icursor != icend; ++icursor) {
- if (colDofMap[feSpaces[j]].find(col(*icursor), colDofIndex) == false)
- continue;
-
- 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 ...
- } // for each row in mat[i][j]
- }
- }
-
-
- if (localMatrix == false) {
- // === 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 - 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)
- nnz.onnz[localRowIdx]++;
- else
- 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++)
- nnz.dnnz[i] = std::min(nnz.dnnz[i], nRankCols);
-
-#if (DEBUG != 0)
- int nMax = 0;
- int nSum = 0;
- for (int i = 0; i < nRankRows; i++) {
- nMax = std::max(nMax, nnz.dnnz[i]);
- nSum += nnz.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
- }
}
diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.h b/AMDiS/src/parallel/PetscSolverGlobalMatrix.h
index 728b1a2e6317ea996539165a27cc7525b6569078..c8624b53c71fc4c678c920f955ad48517a1a4aa6 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.h
+++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.h
@@ -25,55 +25,13 @@
#include <boost/tuple/tuple.hpp>
#include "AMDiS_fwd.h"
+#include "parallel/MatrixNnzStructure.h"
#include "parallel/PetscSolver.h"
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
{
@@ -118,14 +76,8 @@ namespace AMDiS {
virtual void exitPreconditioner(PC pc);
/// Creates a new non zero pattern structure for the PETSc matrix.
- void createPetscNnzStructure(Matrix<DOFMatrix*> *mat);
-
- void createPetscNnzStructureWithCoarseSpace(Matrix<DOFMatrix*> *mat,
- ParallelDofMapping &rowDofMap,
- ParallelDofMapping &colDofMap,
- NnzStructure &nnz,
- bool localMatrix = false);
-
+ // 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);
@@ -148,11 +100,9 @@ namespace AMDiS {
protected:
- NnzStructure nnzInterior;
-
- NnzStructure nnzCoarse;
+ MatrixNnzStructure nnzInterior, nnzCoarse;
- NnzStructure nnzIntCoarse, nnzCoarseInt;
+ MatrixNnzStructure nnzIntCoarse, nnzCoarseInt;
Vec petscSolVec;