42 #ifndef TPETRA_BLOCKCRSMATRIX_HELPERS_DEF_HPP
43 #define TPETRA_BLOCKCRSMATRIX_HELPERS_DEF_HPP
47 #include "Tpetra_BlockCrsMatrix.hpp"
48 #include "Tpetra_CrsMatrix.hpp"
49 #include "Tpetra_HashTable.hpp"
50 #include "Tpetra_Import.hpp"
51 #include "Tpetra_Map.hpp"
52 #include "Tpetra_MultiVector.hpp"
53 #include "Teuchos_ParameterList.hpp"
54 #include "Teuchos_ScalarTraits.hpp"
60 template<
class Scalar,
class LO,
class GO,
class Node>
62 Teuchos::ParameterList pl;
64 out.open(fileName.c_str());
68 template<
class Scalar,
class LO,
class GO,
class Node>
71 out.open(fileName.c_str());
75 template<
class Scalar,
class LO,
class GO,
class Node>
77 Teuchos::ParameterList pl;
81 template<
class Scalar,
class LO,
class GO,
class Node>
87 typedef Teuchos::OrdinalTraits<GO> TOT;
94 RCP<const map_type>
const rowMap = A.
getRowMap();
95 RCP<const Teuchos::Comm<int> > comm = rowMap->getComm();
96 const int myRank = comm->getRank();
97 const size_t numProcs = comm->getSize();
100 bool alwaysUseParallelAlgorithm =
false;
101 if (params.isParameter(
"always use parallel algorithm"))
102 alwaysUseParallelAlgorithm = params.get<
bool>(
"always use parallel algorithm");
103 bool printMatrixMarketHeader =
true;
104 if (params.isParameter(
"print MatrixMarket header"))
105 printMatrixMarketHeader = params.get<
bool>(
"print MatrixMarket header");
107 if (printMatrixMarketHeader && myRank==0) {
108 std::time_t now = std::time(NULL);
110 const std::string dataTypeStr =
111 Teuchos::ScalarTraits<Scalar>::isComplex ?
"complex" :
"real";
121 os <<
"%%MatrixMarket matrix coordinate " << dataTypeStr <<
" general" << std::endl;
122 os <<
"% time stamp: " << ctime(&now);
123 os <<
"% written from " << numProcs <<
" processes" << std::endl;
124 os <<
"% point representation of Tpetra::BlockCrsMatrix" << std::endl;
127 os <<
"% " << numRows <<
" block rows, " << numCols <<
" block columns" << std::endl;
129 os <<
"% block size " << blockSize << std::endl;
130 os << numRows*blockSize <<
" " << numCols*blockSize <<
" " << A.
getGlobalNumEntries()*blockSize*blockSize << std::endl;
133 if (numProcs==1 && !alwaysUseParallelAlgorithm) {
136 size_t numRows = rowMap->getNodeNumElements();
139 RCP<const map_type> allMeshGidsMap = rcp(
new map_type(TOT::invalid(), numRows, A.
getIndexBase(), comm));
142 mv_type allMeshGids(allMeshGidsMap,1);
143 Teuchos::ArrayRCP<GO> allMeshGidsData = allMeshGids.getDataNonConst(0);
145 for (
size_t i=0; i<numRows; i++)
146 allMeshGidsData[i] = rowMap->getGlobalElement(i);
147 allMeshGidsData = Teuchos::null;
150 size_t stripSize = allMeshGids.getGlobalLength() / numProcs;
151 size_t remainder = allMeshGids.getGlobalLength() % numProcs;
153 size_t curStripSize = 0;
154 Teuchos::Array<GO> importMeshGidList;
155 for (
size_t i=0; i<numProcs; i++) {
157 curStripSize = stripSize;
158 if (i<remainder) curStripSize++;
159 importMeshGidList.resize(curStripSize);
160 for (
size_t j=0; j<curStripSize; j++) importMeshGidList[j] = j + curStart + A.
getIndexBase();
161 curStart += curStripSize;
164 TEUCHOS_TEST_FOR_EXCEPTION(myRank>0 && curStripSize!=0,
165 std::runtime_error,
"Tpetra::blockCrsMatrixWriter: (pid "
166 << myRank <<
") map size should be zero, but is " << curStripSize);
167 RCP<map_type> importMeshGidMap = rcp(
new map_type(TOT::invalid(), importMeshGidList(), A.
getIndexBase(), comm));
168 import_type gidImporter(allMeshGidsMap, importMeshGidMap);
169 mv_type importMeshGids(importMeshGidMap, 1);
170 importMeshGids.doImport(allMeshGids, gidImporter,
INSERT);
176 Teuchos::ArrayRCP<const GO> importMeshGidsData = importMeshGids.getData(0);
177 Teuchos::Array<GO> importMeshGidsGO;
178 importMeshGidsGO.reserve(importMeshGidsData.size());
179 for (
typename Teuchos::ArrayRCP<const GO>::size_type j=0; j<importMeshGidsData.size(); ++j)
180 importMeshGidsGO.push_back(importMeshGidsData[j]);
181 RCP<const map_type> importMap = rcp(
new map_type(TOT::invalid(), importMeshGidsGO(), rowMap->getIndexBase(), comm) );
183 import_type importer(rowMap,importMap );
185 RCP<crs_graph_type> graph =
createCrsGraph(importMap,numEntriesPerRow);
186 RCP<const map_type> domainMap = A.getCrsGraph().
getDomainMap();
187 graph->doImport(A.getCrsGraph(), importer,
INSERT);
188 graph->fillComplete(domainMap, importMap);
190 block_crs_matrix_type importA(*graph, A.
getBlockSize());
191 importA.doImport(A, importer,
INSERT);
199 template<
class Scalar,
class LO,
class GO,
class Node>
205 RCP<const map_type> rowMap = A.
getRowMap();
206 RCP<const map_type> colMap = A.
getColMap();
207 RCP<const Teuchos::Comm<int> > comm = rowMap->getComm();
208 const int myRank = comm->getRank();
210 const size_t meshRowOffset = rowMap->getIndexBase();
211 const size_t meshColOffset = colMap->getIndexBase();
212 TEUCHOS_TEST_FOR_EXCEPTION(meshRowOffset != meshColOffset,
213 std::runtime_error,
"Tpetra::writeMatrixStrip: "
214 "mesh row index base != mesh column index base");
219 std::runtime_error,
"Tpetra::writeMatrixStrip: pid "
220 << myRank <<
" should have 0 rows but has " << A.
getNodeNumRows());
222 std::runtime_error,
"Tpetra::writeMatrixStrip: pid "
223 << myRank <<
" should have 0 columns but has " << A.
getNodeNumCols());
228 std::runtime_error,
"Tpetra::writeMatrixStrip: "
229 "number of rows on pid 0 does not match global number of rows");
235 bool precisionChanged=
false;
236 int oldPrecision = 0;
237 if (params.isParameter(
"precision")) {
238 oldPrecision = os.precision(params.get<
int>(
"precision"));
239 precisionChanged=
true;
242 if (params.isParameter(
"zero-based indexing")) {
243 if (params.get<
bool>(
"zero-based indexing") ==
true)
248 for (localRowInd = 0; localRowInd < numLocalRows; ++localRowInd) {
251 const LO* localColInds;
257 GO globalMeshRowID = rowMap->getGlobalElement(localRowInd) - meshRowOffset;
259 for (LO k = 0; k < numEntries; ++k) {
260 GO globalMeshColID = colMap->getGlobalElement(localColInds[k]) - meshColOffset;
261 Scalar*
const curBlock = vals + blockSize * blockSize * k;
263 for (LO j = 0; j < blockSize; ++j) {
264 GO globalPointRowID = globalMeshRowID * blockSize + j + pointOffset;
265 for (LO i = 0; i < blockSize; ++i) {
266 GO globalPointColID = globalMeshColID * blockSize + i + pointOffset;
267 const Scalar curVal = curBlock[i + j * blockSize];
269 os << globalPointRowID <<
" " << globalPointColID <<
" ";
270 if (Teuchos::ScalarTraits<Scalar>::isComplex) {
273 os << Teuchos::ScalarTraits<Scalar>::real (curVal) <<
" "
274 << Teuchos::ScalarTraits<Scalar>::imag (curVal);
284 if (precisionChanged)
285 os.precision(oldPrecision);
286 TEUCHOS_TEST_FOR_EXCEPTION(err != 0,
287 std::runtime_error,
"Tpetra::writeMatrixStrip: "
288 "error getting view of local row " << localRowInd);
294 template<
class Scalar,
class LO,
class GO,
class Node>
295 Teuchos::RCP<BlockCrsMatrix<Scalar, LO, GO, Node> >
307 using Teuchos::Array;
308 using Teuchos::ArrayView;
315 const map_type &pointRowMap = *(pointMatrix.
getRowMap());
316 RCP<const map_type> meshRowMap = createMeshMap<LO,GO,Node>(blockSize, pointRowMap);
318 const map_type &pointColMap = *(pointMatrix.
getColMap());
319 RCP<const map_type> meshColMap = createMeshMap<LO,GO,Node>(blockSize, pointColMap);
321 const map_type &pointDomainMap = *(pointMatrix.
getDomainMap());
322 RCP<const map_type> meshDomainMap = createMeshMap<LO,GO,Node>(blockSize, pointDomainMap);
324 const map_type &pointRangeMap = *(pointMatrix.
getRangeMap());
325 RCP<const map_type> meshRangeMap = createMeshMap<LO,GO,Node>(blockSize, pointRangeMap);
330 RCP<crs_graph_type> meshCrsGraph = rcp(
new crs_graph_type(meshRowMap, meshColMap,
336 ArrayView<const LO> pointColInds;
337 ArrayView<const Scalar> pointVals;
338 Array<GO> meshColGids;
343 for (
int j=0; j<blockSize; ++j) {
344 LO rowLid = i*blockSize+j;
347 for (
int k=0; k<pointColInds.size(); ++k) {
348 GO meshColInd = pointColMap.getGlobalElement(pointColInds[k]) / blockSize;
349 meshColGids.push_back(meshColInd);
354 std::sort(meshColGids.begin(), meshColGids.end());
355 meshColGids.erase( std::unique(meshColGids.begin(), meshColGids.end()), meshColGids.end() );
356 meshCrsGraph->insertGlobalIndices(meshRowMap->getGlobalElement(i), meshColGids());
359 meshCrsGraph->fillComplete(meshDomainMap,meshRangeMap);
362 RCP<block_crs_matrix_type> blockMatrix = rcp(
new block_crs_matrix_type(*meshCrsGraph, blockSize));
365 int maxBlockEntries = blockMatrix->getNodeMaxNumRowEntries();
366 Array<Array<Scalar>> blocks(maxBlockEntries);
367 for (
int i=0; i<maxBlockEntries; ++i)
368 blocks[i].reserve(blockSize*blockSize);
369 std::map<int,int> bcol2bentry;
370 std::map<int,int>::iterator iter;
380 for (
int j=0; j<blockSize; ++j) {
381 LO rowLid = i*blockSize+j;
383 for (
int k=0; k<pointColInds.size(); ++k) {
385 LO meshColInd = pointColInds[k] / blockSize;
386 iter = bcol2bentry.find(meshColInd);
387 if (iter == bcol2bentry.end()) {
389 bcol2bentry[meshColInd] = blkCnt;
390 blocks[blkCnt].push_back(pointVals[k]);
394 int littleBlock = iter->second;
395 blocks[littleBlock].push_back(pointVals[k]);
401 for (iter=bcol2bentry.begin(); iter != bcol2bentry.end(); ++iter) {
402 LO localBlockCol = iter->first;
403 Scalar *vals = (blocks[iter->second]).getRawPtr();
404 blockMatrix->replaceLocalValues(i, &localBlockCol, vals, 1);
408 for (
int j=0; j<maxBlockEntries; ++j)
418 template<
class LO,
class GO,
class Node>
419 Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
422 typedef Teuchos::OrdinalTraits<Tpetra::global_size_t> TOT;
427 Teuchos::Array<GO> meshGids;
433 meshGids.reserve(pointGids.size());
435 for (
int i=0; i<pointGids.size(); ++i) {
436 GO meshGid = (pointGids[i]-indexBase) / blockSize + indexBase;
437 if (hashTable.get(meshGid) == -1) {
438 hashTable.
add(meshGid,1);
439 meshGids.push_back(meshGid);
443 Teuchos::RCP<const map_type> meshMap = Teuchos::rcp(
new map_type(TOT::invalid(), meshGids(), 0, pointMap.
getComm()) );
456 #define TPETRA_BLOCKCRSMATRIX_HELPERS_INSTANT(S,LO,GO,NODE) \
457 template void blockCrsMatrixWriter(BlockCrsMatrix<S,LO,GO,NODE> const &A, std::string const &fileName); \
458 template void blockCrsMatrixWriter(BlockCrsMatrix<S,LO,GO,NODE> const &A, std::string const &fileName, Teuchos::ParameterList const ¶ms); \
459 template void blockCrsMatrixWriter(BlockCrsMatrix<S,LO,GO,NODE> const &A, std::ostream &os); \
460 template void blockCrsMatrixWriter(BlockCrsMatrix<S,LO,GO,NODE> const &A, std::ostream &os, Teuchos::ParameterList const ¶ms); \
461 template void writeMatrixStrip(BlockCrsMatrix<S,LO,GO,NODE> const &A, std::ostream &os, Teuchos::ParameterList const ¶ms); \
462 template Teuchos::RCP<BlockCrsMatrix<S, LO, GO, NODE> > convertToBlockCrsMatrix(const CrsMatrix<S, LO, GO, NODE>& pointMatrix, const LO &blockSize);
467 #define TPETRA_CREATEMESHMAP_INSTANT(LO,GO,NODE) \
468 template Teuchos::RCP<const Map<LO,GO,NODE> > createMeshMap (const LO& blockSize, const Map<LO,GO,NODE>& pointMap);
470 #endif // TPETRA_BLOCKCRSMATRIX_HELPERS_DEF_HPP
Communication plan for data redistribution from a uniquely-owned to a (possibly) multiply-owned distr...
void blockCrsMatrixWriter(BlockCrsMatrix< Scalar, LO, GO, Node > const &A, std::string const &fileName)
Helper function to write a BlockCrsMatrix. Calls the 3-argument version.
Sparse matrix that presents a row-oriented interface that lets users read or modify entries...
Sparse matrix whose entries are small dense square blocks, all of the same dimensions.
size_t getNodeNumRows() const
get the local number of block rows
Teuchos::RCP< const map_type > getRangeMap() const override
The range Map of this matrix.
virtual global_size_t getGlobalNumCols() const
The global number of columns of this matrix.
One or more distributed dense vectors.
Teuchos::RCP< const map_type > getDomainMap() const override
Returns the Map associated with the domain of this graph.
size_t getGlobalMaxNumRowEntries() const override
Maximum number of entries in any row of the graph, over all processes in the graph's communicator...
size_t getNodeNumRows() const override
The number of matrix rows owned by the calling process.
GlobalOrdinal getIndexBase() const
The index base for this Map.
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Accessors for the Teuchos::Comm and Kokkos Node objects.
Teuchos::RCP< CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > createCrsGraph(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &map, size_t maxNumEntriesPerRow=0, const Teuchos::RCP< Teuchos::ParameterList > ¶ms=Teuchos::null)
Nonmember function to create an empty CrsGraph given a row Map and the max number of entries allowed ...
size_t getGlobalMaxNumRowEntries() const override
Maximum number of entries in any row of the matrix, over all processes in the matrix's communicator...
void writeMatrixStrip(BlockCrsMatrix< Scalar, LO, GO, Node > const &A, std::ostream &os, Teuchos::ParameterList const ¶ms)
Helper function called by blockCrsMatrixWriter.
Teuchos::RCP< const Tpetra::Map< LO, GO, Node > > createMeshMap(LO const &blockSize, const Tpetra::Map< LO, GO, Node > &pointMap)
Helper function to generate a mesh map from a point map. Important! It's assumed that point GIDs asso...
virtual size_t getNodeNumCols() const
The number of columns needed to apply the forward operator on this node.
Insert new values that don't currently exist.
Teuchos::RCP< const map_type > getRowMap() const
get the (mesh) map for the rows of this block matrix.
virtual global_size_t getGlobalNumEntries() const
The global number of stored (structurally nonzero) entries.
global_size_t getGlobalNumRows() const
get the global number of block rows
Teuchos::ArrayView< const GlobalOrdinal > getNodeElementList() const
Return a NONOWNING view of the global indices owned by this process.
A distributed graph accessed by rows (adjacency lists) and stored sparsely.
LO getLocalRowView(const LO localRowInd, const LO *&colInds, Scalar *&vals, LO &numInds) const
Get a view of the (mesh, i.e., block) row, using local (mesh, i.e., block) indices.
Teuchos::RCP< const map_type > getDomainMap() const override
The domain Map of this matrix.
void add(const KeyType key, const ValueType value)
Add a key and its value to the hash table.
void getLocalRowView(LocalOrdinal LocalRow, Teuchos::ArrayView< const LocalOrdinal > &indices, Teuchos::ArrayView< const Scalar > &values) const override
Get a constant, nonpersisting view of a row of this matrix, using local row and column indices...
LO getBlockSize() const
The number of degrees of freedom per mesh point.
Teuchos::RCP< BlockCrsMatrix< Scalar, LO, GO, Node > > convertToBlockCrsMatrix(const Tpetra::CrsMatrix< Scalar, LO, GO, Node > &pointMatrix, const LO &blockSize)
Non-member constructor that creates a BlockCrsMatrix from an existing point CrsMatrix.
Teuchos::RCP< const map_type > getColMap() const override
The Map that describes the column distribution in this matrix.
Teuchos::RCP< const map_type > getColMap() const
get the (mesh) map for the columns of this block matrix.
Teuchos::RCP< const map_type > getRowMap() const override
The Map that describes the row distribution in this matrix.
virtual GO getIndexBase() const
The index base for global indices in this matrix.