42 #ifndef TPETRA_EXPERIMENTAL_BLOCKCRSMATRIX_HELPERS_DEF_HPP 43 #define TPETRA_EXPERIMENTAL_BLOCKCRSMATRIX_HELPERS_DEF_HPP 48 #include <Tpetra_HashTable.hpp> 52 namespace Experimental {
54 template<
class Scalar,
class LO,
class GO,
class Node>
56 Teuchos::ParameterList pl;
58 out.open(fileName.c_str());
62 template<
class Scalar,
class LO,
class GO,
class Node>
65 out.open(fileName.c_str());
69 template<
class Scalar,
class LO,
class GO,
class Node>
71 Teuchos::ParameterList pl;
75 template<
class Scalar,
class LO,
class GO,
class Node>
81 typedef Teuchos::OrdinalTraits<GO> TOT;
88 RCP<const map_type>
const rowMap = A.
getRowMap();
89 RCP<const Teuchos::Comm<int> > comm = rowMap->getComm();
90 const int myRank = comm->getRank();
91 const size_t numProcs = comm->getSize();
94 bool alwaysUseParallelAlgorithm =
false;
95 if (params.isParameter(
"always use parallel algorithm"))
96 alwaysUseParallelAlgorithm = params.get<
bool>(
"always use parallel algorithm");
97 bool printMatrixMarketHeader =
true;
98 if (params.isParameter(
"print MatrixMarket header"))
99 printMatrixMarketHeader = params.get<
bool>(
"print MatrixMarket header");
101 if (printMatrixMarketHeader && myRank==0) {
102 std::time_t now = std::time(NULL);
104 const std::string dataTypeStr =
105 Teuchos::ScalarTraits<Scalar>::isComplex ?
"complex" :
"real";
115 os <<
"%%MatrixMarket matrix coordinate " << dataTypeStr <<
" general" << std::endl;
116 os <<
"% time stamp: " << ctime(&now);
117 os <<
"% written from " << numProcs <<
" processes" << std::endl;
118 os <<
"% point representation of Tpetra::Experimental::BlockCrsMatrix" << std::endl;
121 os <<
"% " << numRows <<
" block rows, " << numCols <<
" block columns" << std::endl;
123 os <<
"% block size " << blockSize << std::endl;
124 os << numRows*blockSize <<
" " << numCols*blockSize <<
" " << A.
getGlobalNumEntries()*blockSize*blockSize << std::endl;
127 if (numProcs==1 && !alwaysUseParallelAlgorithm) {
130 size_t numRows = rowMap->getNodeNumElements();
133 RCP<const map_type> allMeshGidsMap = rcp(
new map_type(TOT::invalid(), numRows, A.
getIndexBase(), comm));
136 mv_type allMeshGids(allMeshGidsMap,1);
137 Teuchos::ArrayRCP<GO> allMeshGidsData = allMeshGids.getDataNonConst(0);
139 for (
size_t i=0; i<numRows; i++)
140 allMeshGidsData[i] = rowMap->getGlobalElement(i);
141 allMeshGidsData = Teuchos::null;
144 size_t stripSize = allMeshGids.getGlobalLength() / numProcs;
145 size_t remainder = allMeshGids.getGlobalLength() % numProcs;
147 size_t curStripSize = 0;
148 Teuchos::Array<GO> importMeshGidList;
149 for (
size_t i=0; i<numProcs; i++) {
151 curStripSize = stripSize;
152 if (i<remainder) curStripSize++;
153 importMeshGidList.resize(curStripSize);
154 for (
size_t j=0; j<curStripSize; j++) importMeshGidList[j] = j + curStart + A.
getIndexBase();
155 curStart += curStripSize;
158 TEUCHOS_TEST_FOR_EXCEPTION(myRank>0 && curStripSize!=0,
159 std::runtime_error,
"Tpetra::Experimental::blockCrsMatrixWriter: (pid " 160 << myRank <<
") map size should be zero, but is " << curStripSize);
161 RCP<map_type> importMeshGidMap = rcp(
new map_type(TOT::invalid(), importMeshGidList(), A.
getIndexBase(), comm));
162 import_type gidImporter(allMeshGidsMap, importMeshGidMap);
163 mv_type importMeshGids(importMeshGidMap, 1);
164 importMeshGids.doImport(allMeshGids, gidImporter,
INSERT);
170 Teuchos::ArrayRCP<const GO> importMeshGidsData = importMeshGids.getData(0);
171 Teuchos::Array<GO> importMeshGidsGO;
172 importMeshGidsGO.reserve(importMeshGidsData.size());
173 for (
typename Teuchos::ArrayRCP<const GO>::size_type j=0; j<importMeshGidsData.size(); ++j)
174 importMeshGidsGO.push_back(importMeshGidsData[j]);
175 RCP<const map_type> importMap = rcp(
new map_type(TOT::invalid(), importMeshGidsGO(), rowMap->getIndexBase(), comm) );
177 import_type importer(rowMap,importMap );
179 RCP<crs_graph_type> graph =
createCrsGraph(importMap,numEntriesPerRow);
180 RCP<const map_type> domainMap = A.getCrsGraph().
getDomainMap();
181 graph->doImport(A.getCrsGraph(), importer,
INSERT);
182 graph->fillComplete(domainMap, importMap);
184 block_crs_matrix_type importA(*graph, A.
getBlockSize());
185 importA.doImport(A, importer,
INSERT);
193 template<
class Scalar,
class LO,
class GO,
class Node>
199 RCP<const map_type> rowMap = A.
getRowMap();
200 RCP<const map_type> colMap = A.
getColMap();
201 RCP<const Teuchos::Comm<int> > comm = rowMap->getComm();
202 const int myRank = comm->getRank();
204 const size_t meshRowOffset = rowMap->getIndexBase();
205 const size_t meshColOffset = colMap->getIndexBase();
206 TEUCHOS_TEST_FOR_EXCEPTION(meshRowOffset != meshColOffset,
207 std::runtime_error,
"Tpetra::Experimental::writeMatrixStrip: " 208 "mesh row index base != mesh column index base");
213 std::runtime_error,
"Tpetra::Experimental::writeMatrixStrip: pid " 214 << myRank <<
" should have 0 rows but has " << A.
getNodeNumRows());
216 std::runtime_error,
"Tpetra::Experimental::writeMatrixStrip: pid " 217 << myRank <<
" should have 0 columns but has " << A.
getNodeNumCols());
222 std::runtime_error,
"Tpetra::Experimental::writeMatrixStrip: " 223 "number of rows on pid 0 does not match global number of rows");
229 bool precisionChanged=
false;
230 int oldPrecision = 0;
231 if (params.isParameter(
"precision")) {
232 oldPrecision = os.precision(params.get<
int>(
"precision"));
233 precisionChanged=
true;
236 if (params.isParameter(
"zero-based indexing")) {
237 if (params.get<
bool>(
"zero-based indexing") ==
true)
242 for (localRowInd = 0; localRowInd < numLocalRows; ++localRowInd) {
245 const LO* localColInds;
251 GO globalMeshRowID = rowMap->getGlobalElement(localRowInd) - meshRowOffset;
253 for (LO k = 0; k < numEntries; ++k) {
254 GO globalMeshColID = colMap->getGlobalElement(localColInds[k]) - meshColOffset;
255 Scalar*
const curBlock = vals + blockSize * blockSize * k;
257 for (LO j = 0; j < blockSize; ++j) {
258 GO globalPointRowID = globalMeshRowID * blockSize + j + pointOffset;
259 for (LO i = 0; i < blockSize; ++i) {
260 GO globalPointColID = globalMeshColID * blockSize + i + pointOffset;
261 const Scalar curVal = curBlock[i + j * blockSize];
263 os << globalPointRowID <<
" " << globalPointColID <<
" ";
264 if (Teuchos::ScalarTraits<Scalar>::isComplex) {
267 os << Teuchos::ScalarTraits<Scalar>::real (curVal) <<
" " 268 << Teuchos::ScalarTraits<Scalar>::imag (curVal);
278 if (precisionChanged)
279 os.precision(oldPrecision);
280 TEUCHOS_TEST_FOR_EXCEPTION(err != 0,
281 std::runtime_error,
"Tpetra::Experimental::writeMatrixStrip: " 282 "error getting view of local row " << localRowInd);
288 template<
class Scalar,
class LO,
class GO,
class Node>
289 Teuchos::RCP<BlockCrsMatrix<Scalar, LO, GO, Node> >
301 using Teuchos::Array;
302 using Teuchos::ArrayView;
309 const map_type &pointRowMap = *(pointMatrix.
getRowMap());
310 RCP<const map_type> meshRowMap = createMeshMap<LO,GO,Node>(blockSize, pointRowMap);
312 const map_type &pointColMap = *(pointMatrix.
getColMap());
313 RCP<const map_type> meshColMap = createMeshMap<LO,GO,Node>(blockSize, pointColMap);
315 const map_type &pointDomainMap = *(pointMatrix.
getDomainMap());
316 RCP<const map_type> meshDomainMap = createMeshMap<LO,GO,Node>(blockSize, pointDomainMap);
318 const map_type &pointRangeMap = *(pointMatrix.
getRangeMap());
319 RCP<const map_type> meshRangeMap = createMeshMap<LO,GO,Node>(blockSize, pointRangeMap);
324 RCP<crs_graph_type> meshCrsGraph = rcp(
new crs_graph_type(meshRowMap, meshColMap,
330 ArrayView<const LO> pointColInds;
331 ArrayView<const Scalar> pointVals;
332 Array<GO> meshColGids;
337 for (
int j=0; j<blockSize; ++j) {
338 LO rowLid = i*blockSize+j;
341 for (
int k=0; k<pointColInds.size(); ++k) {
342 GO meshColInd = pointColMap.getGlobalElement(pointColInds[k]) / blockSize;
343 meshColGids.push_back(meshColInd);
348 std::sort(meshColGids.begin(), meshColGids.end());
349 meshColGids.erase( std::unique(meshColGids.begin(), meshColGids.end()), meshColGids.end() );
350 meshCrsGraph->insertGlobalIndices(meshRowMap->getGlobalElement(i), meshColGids());
353 meshCrsGraph->fillComplete(meshDomainMap,meshRangeMap);
356 RCP<block_crs_matrix_type> blockMatrix = rcp(
new block_crs_matrix_type(*meshCrsGraph, blockSize));
359 int maxBlockEntries = blockMatrix->getNodeMaxNumRowEntries();
360 Array<Array<Scalar>> blocks(maxBlockEntries);
361 for (
int i=0; i<maxBlockEntries; ++i)
362 blocks[i].reserve(blockSize*blockSize);
363 std::map<int,int> bcol2bentry;
364 std::map<int,int>::iterator iter;
374 for (
int j=0; j<blockSize; ++j) {
375 LO rowLid = i*blockSize+j;
377 for (
int k=0; k<pointColInds.size(); ++k) {
379 LO meshColInd = pointColInds[k] / blockSize;
380 iter = bcol2bentry.find(meshColInd);
381 if (iter == bcol2bentry.end()) {
383 bcol2bentry[meshColInd] = blkCnt;
384 blocks[blkCnt].push_back(pointVals[k]);
388 int littleBlock = iter->second;
389 blocks[littleBlock].push_back(pointVals[k]);
395 for (iter=bcol2bentry.begin(); iter != bcol2bentry.end(); ++iter) {
396 LO localBlockCol = iter->first;
397 Scalar *vals = (blocks[iter->second]).getRawPtr();
398 blockMatrix->replaceLocalValues(i, &localBlockCol, vals, 1);
402 for (
int j=0; j<maxBlockEntries; ++j)
412 template<
class LO,
class GO,
class Node>
413 Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
416 typedef Teuchos::OrdinalTraits<Tpetra::global_size_t> TOT;
421 Teuchos::Array<GO> meshGids;
427 meshGids.reserve(pointGids.size());
429 for (
int i=0; i<pointGids.size(); ++i) {
430 GO meshGid = (pointGids[i]-indexBase) / blockSize + indexBase;
431 if (hashTable.get(meshGid) == -1) {
432 hashTable.
add(meshGid,1);
433 meshGids.push_back(meshGid);
437 Teuchos::RCP<const map_type> meshMap = Teuchos::rcp(
new map_type(TOT::invalid(), meshGids(), 0, pointMap.
getComm()) );
451 #define TPETRA_EXPERIMENTAL_BLOCKCRSMATRIX_HELPERS_INSTANT(S,LO,GO,NODE) \ 452 template void Experimental::blockCrsMatrixWriter(Experimental::BlockCrsMatrix<S,LO,GO,NODE> const &A, std::string const &fileName); \ 453 template void Experimental::blockCrsMatrixWriter(Experimental::BlockCrsMatrix<S,LO,GO,NODE> const &A, std::string const &fileName, Teuchos::ParameterList const ¶ms); \ 454 template void Experimental::blockCrsMatrixWriter(Experimental::BlockCrsMatrix<S,LO,GO,NODE> const &A, std::ostream &os); \ 455 template void Experimental::blockCrsMatrixWriter(Experimental::BlockCrsMatrix<S,LO,GO,NODE> const &A, std::ostream &os, Teuchos::ParameterList const ¶ms); \ 456 template void Experimental::writeMatrixStrip(Experimental::BlockCrsMatrix<S,LO,GO,NODE> const &A, std::ostream &os, Teuchos::ParameterList const ¶ms); \ 457 template Teuchos::RCP<Experimental::BlockCrsMatrix<S, LO, GO, NODE> > Experimental::convertToBlockCrsMatrix(const CrsMatrix<S, LO, GO, NODE>& pointMatrix, const LO &blockSize); 464 #define TPETRA_EXPERIMENTAL_CREATEMESHMAP_INSTANT(LO,GO,NODE) \ 465 template Teuchos::RCP<const Map<LO,GO,NODE> > createMeshMap (const LO& blockSize, const Map<LO,GO,NODE>& pointMap); 467 #endif // TPETRA_EXPERIMENTAL_BLOCKCRSMATRIX_HELPERS_DEF_HPP 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.
Communication plan for data redistribution from a uniquely-owned to a (possibly) multiply-owned distr...
Namespace Tpetra contains the class and methods constituting the Tpetra library.
Sparse matrix that presents a row-oriented interface that lets users read or modify entries...
global_size_t getGlobalNumRows() const
get the global number of block rows
void writeMatrixStrip(BlockCrsMatrix< Scalar, LO, GO, Node > const &A, std::ostream &os, Teuchos::ParameterList const ¶ms)
Helper function called by blockCrsMatrixWriter.
Teuchos::RCP< const map_type > getRowMap() const
The Map that describes the row distribution in this matrix.
One or more distributed dense vectors.
virtual global_size_t getGlobalNumCols() const
The global number of columns of this matrix.
Teuchos::RCP< CrsGraph< LocalOrdinal, GlobalOrdinal, Node, classic > > 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 ...
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.
void blockCrsMatrixWriter(BlockCrsMatrix< Scalar, LO, GO, Node > const &A, std::string const &fileName)
Helper function to write a BlockCrsMatrix. Calls the 3-argument version.
Teuchos::RCP< const map_type > getDomainMap() const
Returns the Map associated with the domain of this graph.
virtual size_t getNodeNumCols() const
The number of columns needed to apply the forward operator on this node.
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Accessors for the Teuchos::Comm and Kokkos Node objects.
Insert new values that don't currently exist.
GlobalOrdinal getIndexBase() const
The index base for this Map.
void getLocalRowView(LocalOrdinal LocalRow, Teuchos::ArrayView< const LocalOrdinal > &indices, Teuchos::ArrayView< const Scalar > &values) const
Get a constant, nonpersisting view of a row of this matrix, using local row and column indices...
Teuchos::RCP< const map_type > getRowMap() const
get the (mesh) map for the rows of this block matrix.
Teuchos::ArrayView< const GlobalOrdinal > getNodeElementList() const
Return a NONOWNING view of the global indices owned by this process.
size_t getGlobalMaxNumRowEntries() const
Maximum number of entries in all rows over all processes.
size_t getNodeNumRows() const
The number of matrix rows owned by the calling process.
A distributed graph accessed by rows (adjacency lists) and stored sparsely.
virtual global_size_t getGlobalNumEntries() const
The global number of stored (structurally nonzero) entries.
void add(const KeyType key, const ValueType value)
Add a key and its value to the hash table.
virtual GO getIndexBase() const
The index base for global indices in this matrix.
size_t getGlobalMaxNumRowEntries() const
Returns the maximum number of entries across all rows/columns on all nodes.
Sparse matrix whose entries are small dense square blocks, all of the same dimensions.
Teuchos::RCP< const map_type > getColMap() const
The Map that describes the column distribution in this matrix.
Teuchos::RCP< const map_type > getDomainMap() const
The domain Map of this matrix.
LO getBlockSize() const
The number of degrees of freedom per mesh point.
Teuchos::RCP< const map_type > getRangeMap() const
The range Map of this matrix.
Teuchos::RCP< const map_type > getColMap() const
get the (mesh) map for the columns of this block matrix.
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...
size_t getNodeNumRows() const
get the local number of block rows