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 RCP<BlockCrsMatrix<Scalar, LO, GO, Node>>
301 using Teuchos::ArrayView;
302 using Teuchos::Array;
308 const map_type &pointRowMap = *(pointMatrix.
getRowMap());
309 RCP<const map_type> meshRowMap = createMeshMap<LO,GO,Node>(blockSize, pointRowMap);
311 const map_type &pointColMap = *(pointMatrix.
getColMap());
312 RCP<const map_type> meshColMap = createMeshMap<LO,GO,Node>(blockSize, pointColMap);
314 const map_type &pointDomainMap = *(pointMatrix.
getDomainMap());
315 RCP<const map_type> meshDomainMap = createMeshMap<LO,GO,Node>(blockSize, pointDomainMap);
317 const map_type &pointRangeMap = *(pointMatrix.
getRangeMap());
318 RCP<const map_type> meshRangeMap = createMeshMap<LO,GO,Node>(blockSize, pointRangeMap);
323 RCP<crs_graph_type> meshCrsGraph = rcp(
new crs_graph_type(meshRowMap, meshColMap,
329 ArrayView<const LO> pointColInds;
330 ArrayView<const Scalar> pointVals;
331 Array<GO> meshColGids;
336 for (
int j=0; j<blockSize; ++j) {
337 LO rowLid = i*blockSize+j;
340 for (
int k=0; k<pointColInds.size(); ++k) {
341 GO meshColInd = pointColMap.getGlobalElement(pointColInds[k]) / blockSize;
342 meshColGids.push_back(meshColInd);
347 std::sort(meshColGids.begin(), meshColGids.end());
348 meshColGids.erase( std::unique(meshColGids.begin(), meshColGids.end()), meshColGids.end() );
349 meshCrsGraph->insertGlobalIndices(meshRowMap->getGlobalElement(i), meshColGids());
352 meshCrsGraph->fillComplete(meshDomainMap,meshRangeMap);
355 RCP<block_crs_matrix_type> blockMatrix = rcp(
new block_crs_matrix_type(*meshCrsGraph, blockSize));
358 int maxBlockEntries = blockMatrix->getNodeMaxNumRowEntries();
359 Array<Array<Scalar>> blocks(maxBlockEntries);
360 for (
int i=0; i<maxBlockEntries; ++i)
361 blocks[i].reserve(blockSize*blockSize);
362 std::map<int,int> bcol2bentry;
363 std::map<int,int>::iterator iter;
373 for (
int j=0; j<blockSize; ++j) {
374 LO rowLid = i*blockSize+j;
376 for (
int k=0; k<pointColInds.size(); ++k) {
378 LO meshColInd = pointColInds[k] / blockSize;
379 iter = bcol2bentry.find(meshColInd);
380 if (iter == bcol2bentry.end()) {
382 bcol2bentry[meshColInd] = blkCnt;
383 blocks[blkCnt].push_back(pointVals[k]);
387 int littleBlock = iter->second;
388 blocks[littleBlock].push_back(pointVals[k]);
394 for (iter=bcol2bentry.begin(); iter != bcol2bentry.end(); ++iter) {
395 LO localBlockCol = iter->first;
396 Scalar *vals = (blocks[iter->second]).getRawPtr();
397 blockMatrix->replaceLocalValues(i, &localBlockCol, vals, 1);
401 for (
int j=0; j<maxBlockEntries; ++j)
411 template<
class LO,
class GO,
class Node>
412 Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
415 typedef Teuchos::OrdinalTraits<Tpetra::global_size_t> TOT;
420 Teuchos::Array<GO> meshGids;
426 meshGids.reserve(pointGids.size());
428 for (
int i=0; i<pointGids.size(); ++i) {
429 GO meshGid = (pointGids[i]-indexBase) / blockSize + indexBase;
430 if (hashTable.get(meshGid) == -1) {
431 hashTable.
add(meshGid,1);
432 meshGids.push_back(meshGid);
436 Teuchos::RCP<const map_type> meshMap = Teuchos::rcp(
new map_type(TOT::invalid(), meshGids(), 0, pointMap.
getComm()) );
450 #define TPETRA_EXPERIMENTAL_BLOCKCRSMATRIX_HELPERS_INSTANT(S,LO,GO,NODE) \ 451 template void Experimental::blockCrsMatrixWriter(Experimental::BlockCrsMatrix<S,LO,GO,NODE> const &A, std::string const &fileName); \ 452 template void Experimental::blockCrsMatrixWriter(Experimental::BlockCrsMatrix<S,LO,GO,NODE> const &A, std::string const &fileName, Teuchos::ParameterList const ¶ms); \ 453 template void Experimental::blockCrsMatrixWriter(Experimental::BlockCrsMatrix<S,LO,GO,NODE> const &A, std::ostream &os); \ 454 template void Experimental::blockCrsMatrixWriter(Experimental::BlockCrsMatrix<S,LO,GO,NODE> const &A, std::ostream &os, Teuchos::ParameterList const ¶ms); \ 455 template void Experimental::writeMatrixStrip(Experimental::BlockCrsMatrix<S,LO,GO,NODE> const &A, std::ostream &os, Teuchos::ParameterList const ¶ms); \ 456 template Teuchos::RCP<Experimental::BlockCrsMatrix<S, LO, GO, NODE> > Experimental::convertToBlockCrsMatrix(const CrsMatrix<S, LO, GO, NODE>& pointMatrix, const LO &blockSize); 463 #define TPETRA_EXPERIMENTAL_CREATEMESHMAP_INSTANT(LO,GO,NODE) \ 464 template Teuchos::RCP<const Map<LO,GO,NODE> > createMeshMap (const LO& blockSize, const Map<LO,GO,NODE>& pointMap); 466 #endif // TPETRA_EXPERIMENTAL_BLOCKCRSMATRIX_HELPERS_DEF_HPP 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...
Teuchos::RCP< const map_type > getColMap() const
get the (mesh) map for the columns of this block matrix.
LO getBlockSize() const
The number of degrees of freedom per mesh point.
Teuchos::RCP< const map_type > getRowMap() const
get the (mesh) map for the rows of this block matrix.
void writeMatrixStrip(BlockCrsMatrix< Scalar, LO, GO, Node > const &A, std::ostream &os, Teuchos::ParameterList const ¶ms)
Helper function called by blockCrsMatrixWriter.
size_t getGlobalMaxNumRowEntries() const
Maximum number of entries in all rows over all processes.
Teuchos::RCP< const map_type > getRowMap() const
The Map that describes the row distribution in this matrix.
size_t getNodeNumRows() const
get the local number of block rows
Teuchos::RCP< const map_type > getColMap() const
The Map that describes the column distribution in this matrix.
One or more distributed dense vectors.
global_size_t getGlobalNumRows() const
get the global number of block rows
GlobalOrdinal getIndexBase() const
The index base for this Map.
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 ...
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Accessors for the Teuchos::Comm and Kokkos Node objects.
void blockCrsMatrixWriter(BlockCrsMatrix< Scalar, LO, GO, Node > const &A, std::string const &fileName)
Helper function to write a BlockCrsMatrix. Calls the 3-argument version.
size_t getNodeNumRows() const
The number of matrix rows owned by the calling process.
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...
virtual GO getIndexBase() const
The index base for global indices in this matrix.
Insert new values that don't currently exist.
Teuchos::RCP< const map_type > getDomainMap() const
Returns the Map associated with the domain of this graph.
virtual global_size_t getGlobalNumCols() const
The global number of columns of this matrix.
Teuchos::RCP< const map_type > getDomainMap() const
The domain Map of this matrix.
virtual size_t getNodeNumCols() const
The number of columns needed to apply the forward operator on this node.
virtual global_size_t getGlobalNumEntries() const
The global number of stored (structurally nonzero) entries.
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...
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::ArrayView< const GlobalOrdinal > getNodeElementList() const
Return a view of the global indices owned by this process.
A distributed graph accessed by rows (adjacency lists) and stored sparsely.
size_t getGlobalMaxNumRowEntries() const
Returns the maximum number of entries across all rows/columns on all nodes.
void add(const KeyType key, const ValueType value)
Add a key and its value to the hash table.
Constant block CRS matrix class.
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 > getRangeMap() const
The range Map of this matrix.