Tpetra parallel linear algebra Version of the Day
Loading...
Searching...
No Matches
Tpetra_BlockCrsMatrix_Helpers_def.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Tpetra: Templated Linear Algebra Services Package
5// Copyright (2008) Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// ************************************************************************
38// @HEADER
39
40#ifndef TPETRA_BLOCKCRSMATRIX_HELPERS_DEF_HPP
41#define TPETRA_BLOCKCRSMATRIX_HELPERS_DEF_HPP
42
44
45#include "Tpetra_BlockCrsMatrix.hpp"
46#include "Tpetra_CrsMatrix.hpp"
47#include "Tpetra_HashTable.hpp"
48#include "Tpetra_Import.hpp"
49#include "Tpetra_Map.hpp"
50#include "Tpetra_MultiVector.hpp"
51#include "Teuchos_ParameterList.hpp"
52#include "Teuchos_ScalarTraits.hpp"
53#include <ctime>
54#include <fstream>
55
56namespace Tpetra {
57
58 template<class Scalar, class LO, class GO, class Node>
60 Teuchos::ParameterList pl;
61 std::ofstream out;
62 out.open(fileName.c_str());
64 }
65
66 template<class Scalar, class LO, class GO, class Node>
67 void blockCrsMatrixWriter(BlockCrsMatrix<Scalar,LO,GO,Node> const &A, std::string const &fileName, Teuchos::ParameterList const &params) {
68 std::ofstream out;
69 out.open(fileName.c_str());
71 }
72
73 template<class Scalar, class LO, class GO, class Node>
75 Teuchos::ParameterList pl;
77 }
78
79 template<class Scalar, class LO, class GO, class Node>
80 void blockCrsMatrixWriter(BlockCrsMatrix<Scalar,LO,GO,Node> const &A, std::ostream &os, Teuchos::ParameterList const &params) {
81
82 using Teuchos::RCP;
83 using Teuchos::rcp;
84
85 typedef Teuchos::OrdinalTraits<GO> TOT;
86 typedef BlockCrsMatrix<Scalar, LO, GO, Node> block_crs_matrix_type;
87 typedef Tpetra::Import<LO, GO, Node> import_type;
88 typedef Tpetra::Map<LO, GO, Node> map_type;
90 typedef Tpetra::CrsGraph<LO, GO, Node> crs_graph_type;
91
92 RCP<const map_type> const rowMap = A.getRowMap(); //"mesh" map
93 RCP<const Teuchos::Comm<int> > comm = rowMap->getComm();
94 const int myRank = comm->getRank();
95 const size_t numProcs = comm->getSize();
96
97 // If true, force use of the import strip-mining infrastructure. This is useful for debugging on one process.
98 bool alwaysUseParallelAlgorithm = false;
99 if (params.isParameter("always use parallel algorithm"))
100 alwaysUseParallelAlgorithm = params.get<bool>("always use parallel algorithm");
101 bool printMatrixMarketHeader = true;
102 if (params.isParameter("print MatrixMarket header"))
103 printMatrixMarketHeader = params.get<bool>("print MatrixMarket header");
104
106 std::time_t now = std::time(NULL);
107
108 const std::string dataTypeStr =
109 Teuchos::ScalarTraits<Scalar>::isComplex ? "complex" : "real";
110
111 // Explanation of first line of file:
112 // - "%%MatrixMarket" is the tag for Matrix Market format.
113 // - "matrix" is what we're printing.
114 // - "coordinate" means sparse (triplet format), rather than dense.
115 // - "real" / "complex" is the type (in an output format sense,
116 // not in a C++ sense) of each value of the matrix. (The
117 // other two possibilities are "integer" (not really necessary
118 // here) and "pattern" (no values, just graph).)
119 os << "%%MatrixMarket matrix coordinate " << dataTypeStr << " general" << std::endl;
120 os << "% time stamp: " << ctime(&now);
121 os << "% written from " << numProcs << " processes" << std::endl;
122 os << "% point representation of Tpetra::BlockCrsMatrix" << std::endl;
123 size_t numRows = A.getGlobalNumRows();
124 size_t numCols = A.getGlobalNumCols();
125 os << "% " << numRows << " block rows, " << numCols << " block columns" << std::endl;
126 const LO blockSize = A.getBlockSize();
127 os << "% block size " << blockSize << std::endl;
128 os << numRows*blockSize << " " << numCols*blockSize << " " << A.getGlobalNumEntries()*blockSize*blockSize << std::endl;
129 }
130
133 } else {
134 size_t numRows = rowMap->getLocalNumElements();
135
136 //Create source map
137 RCP<const map_type> allMeshGidsMap = rcp(new map_type(TOT::invalid(), numRows, A.getIndexBase(), comm));
138 //Create and populate vector of mesh GIDs corresponding to this pid's rows.
139 //This vector will be imported one pid's worth of information at a time to pid 0.
140 mv_type allMeshGids(allMeshGidsMap,1);
141 Teuchos::ArrayRCP<GO> allMeshGidsData = allMeshGids.getDataNonConst(0);
142
143 for (size_t i=0; i<numRows; i++)
144 allMeshGidsData[i] = rowMap->getGlobalElement(i);
145 allMeshGidsData = Teuchos::null;
146
147 // Now construct a RowMatrix on PE 0 by strip-mining the rows of the input matrix A.
148 size_t stripSize = allMeshGids.getGlobalLength() / numProcs;
149 size_t remainder = allMeshGids.getGlobalLength() % numProcs;
150 size_t curStart = 0;
151 size_t curStripSize = 0;
152 Teuchos::Array<GO> importMeshGidList;
153 for (size_t i=0; i<numProcs; i++) {
154 if (myRank==0) { // Only PE 0 does this part
156 if (i<remainder) curStripSize++; // handle leftovers
157 importMeshGidList.resize(curStripSize); // Set size of vector to max needed
158 for (size_t j=0; j<curStripSize; j++) importMeshGidList[j] = j + curStart + A.getIndexBase();
160 }
161 // The following import map should be non-trivial only on PE 0.
163 std::runtime_error, "Tpetra::blockCrsMatrixWriter: (pid "
164 << myRank << ") map size should be zero, but is " << curStripSize);
165 RCP<map_type> importMeshGidMap = rcp(new map_type(TOT::invalid(), importMeshGidList(), A.getIndexBase(), comm));
169
170 // importMeshGids now has a list of GIDs for the current strip of matrix rows.
171 // Use these values to build another importer that will get rows of the matrix.
172
173 // The following import map will be non-trivial only on PE 0.
174 Teuchos::ArrayRCP<const GO> importMeshGidsData = importMeshGids.getData(0);
175 Teuchos::Array<GO> importMeshGidsGO;
176 importMeshGidsGO.reserve(importMeshGidsData.size());
177 for (typename Teuchos::ArrayRCP<const GO>::size_type j=0; j<importMeshGidsData.size(); ++j)
179 RCP<const map_type> importMap = rcp(new map_type(TOT::invalid(), importMeshGidsGO(), rowMap->getIndexBase(), comm) );
180
181 import_type importer(rowMap,importMap );
182 size_t numEntriesPerRow = A.getCrsGraph().getGlobalMaxNumRowEntries();
184 RCP<const map_type> domainMap = A.getCrsGraph().getDomainMap();
185 graph->doImport(A.getCrsGraph(), importer, INSERT);
186 graph->fillComplete(domainMap, importMap);
187
188 block_crs_matrix_type importA(*graph, A.getBlockSize());
189 importA.doImport(A, importer, INSERT);
190
191 // Finally we are ready to write this strip of the matrix
193 }
194 }
195 }
196
197 template<class Scalar, class LO, class GO, class Node>
198 void writeMatrixStrip(BlockCrsMatrix<Scalar,LO,GO,Node> const &A, std::ostream &os, Teuchos::ParameterList const &params) {
199 using Teuchos::RCP;
200 using map_type = Tpetra::Map<LO, GO, Node>;
202 using bcrs_local_inds_host_view_type = typename bcrs_type::local_inds_host_view_type;
203 using bcrs_values_host_view_type = typename bcrs_type::values_host_view_type;
204 using impl_scalar_type = typename bcrs_type::impl_scalar_type;
205
206 size_t numRows = A.getGlobalNumRows();
207 RCP<const map_type> rowMap = A.getRowMap();
208 RCP<const map_type> colMap = A.getColMap();
209 RCP<const Teuchos::Comm<int> > comm = rowMap->getComm();
210 const int myRank = comm->getRank();
211
212 const size_t meshRowOffset = rowMap->getIndexBase();
213 const size_t meshColOffset = colMap->getIndexBase();
215 std::runtime_error, "Tpetra::writeMatrixStrip: "
216 "mesh row index base != mesh column index base");
217
218 if (myRank !=0) {
219
220 TEUCHOS_TEST_FOR_EXCEPTION(A.getLocalNumRows() != 0,
221 std::runtime_error, "Tpetra::writeMatrixStrip: pid "
222 << myRank << " should have 0 rows but has " << A.getLocalNumRows());
223 TEUCHOS_TEST_FOR_EXCEPTION(A.getLocalNumCols() != 0,
224 std::runtime_error, "Tpetra::writeMatrixStrip: pid "
225 << myRank << " should have 0 columns but has " << A.getLocalNumCols());
226
227 } else {
228
229 TEUCHOS_TEST_FOR_EXCEPTION(numRows != A.getLocalNumRows(),
230 std::runtime_error, "Tpetra::writeMatrixStrip: "
231 "number of rows on pid 0 does not match global number of rows");
232
233
234 int err = 0;
235 const LO blockSize = A.getBlockSize();
236 const size_t numLocalRows = A.getLocalNumRows();
237 bool precisionChanged=false;
238 int oldPrecision = 0; // avoid "unused variable" warning
239 if (params.isParameter("precision")) {
240 oldPrecision = os.precision(params.get<int>("precision"));
241 precisionChanged=true;
242 }
243 int pointOffset = 1;
244 if (params.isParameter("zero-based indexing")) {
245 if (params.get<bool>("zero-based indexing") == true)
246 pointOffset = 0;
247 }
248
249 size_t localRowInd;
251
252 // Get a view of the current row.
255 LO numEntries;
256 A.getLocalRowView (localRowInd, localColInds, vals); numEntries = localColInds.extent(0);
257 GO globalMeshRowID = rowMap->getGlobalElement(localRowInd) - meshRowOffset;
258
259 for (LO k = 0; k < numEntries; ++k) {
260 GO globalMeshColID = colMap->getGlobalElement(localColInds[k]) - meshColOffset;
261 const impl_scalar_type* curBlock = vals.data() + blockSize * blockSize * k;
262 // Blocks are stored in row-major format.
263 for (LO j = 0; j < blockSize; ++j) {
265 for (LO i = 0; i < blockSize; ++i) {
267 const impl_scalar_type curVal = curBlock[i + j * blockSize];
268
269 os << globalPointRowID << " " << globalPointColID << " ";
270 if (Teuchos::ScalarTraits<impl_scalar_type>::isComplex) {
271 // Matrix Market format wants complex values to be
272 // written as space-delimited pairs. See Bug 6469.
274 << Teuchos::ScalarTraits<impl_scalar_type>::imag (curVal);
275 }
276 else {
277 os << curVal;
278 }
279 os << std::endl;
280 }
281 }
282 }
283 }
285 os.precision(oldPrecision);
287 std::runtime_error, "Tpetra::writeMatrixStrip: "
288 "error getting view of local row " << localRowInd);
289
290 }
291
292 }
293
294 template<class Scalar, class LO, class GO, class Node>
295 Teuchos::RCP<BlockCrsMatrix<Scalar, LO, GO, Node> >
297 {
298
299 /*
300 ASSUMPTIONS:
301
302 1) In point matrix, all entries associated with a little block are present (even if they are zero).
303 2) For given mesh DOF, point DOFs appear consecutively and in ascending order in row & column maps.
304 3) Point column map and block column map are ordered consistently.
305 */
306
307 using Teuchos::Array;
308 using Teuchos::ArrayView;
309 using Teuchos::RCP;
310
311 typedef Tpetra::BlockCrsMatrix<Scalar,LO,GO,Node> block_crs_matrix_type;
312 typedef Tpetra::Map<LO,GO,Node> map_type;
313 typedef Tpetra::CrsGraph<LO,GO,Node> crs_graph_type;
314 typedef Tpetra::CrsMatrix<Scalar, LO,GO,Node> crs_matrix_type;
315
316 const map_type &pointRowMap = *(pointMatrix.getRowMap());
318
319 const map_type &pointColMap = *(pointMatrix.getColMap());
321 if(meshColMap.is_null()) throw std::runtime_error("ERROR: Cannot create mesh colmap");
322
323 const map_type &pointDomainMap = *(pointMatrix.getDomainMap());
325
326 const map_type &pointRangeMap = *(pointMatrix.getRangeMap());
328
329 // Use graph ctor that provides column map and upper bound on nonzeros per row.
330 // We can use static profile because the point graph should have at least as many entries per
331 // row as the mesh graph.
333 pointMatrix.getGlobalMaxNumRowEntries()));
334 // Fill the graph by walking through the matrix. For each mesh row, we query the collection of point
335 // rows associated with it. The point column ids are converted to mesh column ids and put into an array.
336 // As each point row collection is finished, the mesh column ids are sorted, made unique, and inserted
337 // into the mesh graph.
338 typename crs_matrix_type::local_inds_host_view_type pointColInds;
339 typename crs_matrix_type::values_host_view_type pointVals;
341 meshColGids.reserve(pointMatrix.getGlobalMaxNumRowEntries());
342
343 //again, I assume that point GIDs associated with a mesh GID are consecutive.
344 //if they are not, this will break!!
345 GO indexBase = pointColMap.getIndexBase();
346 for (size_t i=0; i<pointMatrix.getLocalNumRows()/blockSize; i++) {
347 for (int j=0; j<blockSize; ++j) {
348 LO rowLid = i*blockSize+j;
349 pointMatrix.getLocalRowView(rowLid,pointColInds,pointVals); //TODO optimization: Since I don't care about values,
350 //TODO I should use the graph instead.
351 for (size_t k=0; k<pointColInds.size(); ++k) {
352 GO meshColInd = (pointColMap.getGlobalElement(pointColInds[k]) - indexBase) / blockSize + indexBase;
353 if (meshColMap->getLocalElement(meshColInd) == Teuchos::OrdinalTraits<GO>::invalid()) {
354 std::ostringstream oss;
355 oss<< "["<<i<<"] ERROR: meshColId "<< meshColInd <<" is not in the column map. Correspnds to pointColId = "<<pointColInds[k]<<std::endl;
356 throw std::runtime_error(oss.str());
357 }
358
359 meshColGids.push_back(meshColInd);
360 }
361 }
362 //List of mesh GIDs probably contains duplicates because we looped over all point rows in the block.
363 //Sort and make unique.
364 std::sort(meshColGids.begin(), meshColGids.end());
365 meshColGids.erase( std::unique(meshColGids.begin(), meshColGids.end()), meshColGids.end() );
366 meshCrsGraph->insertGlobalIndices(meshRowMap->getGlobalElement(i), meshColGids());
367 meshColGids.clear();
368 }
370
371 //create and populate the block matrix
372 RCP<block_crs_matrix_type> blockMatrix = rcp(new block_crs_matrix_type(*meshCrsGraph, blockSize));
373
376
377 //preallocate the maximum number of (dense) block entries needed by any row
378 int maxBlockEntries = blockMatrix->getLocalMaxNumRowEntries();
380 for (int i=0; i<maxBlockEntries; ++i)
382 std::map<int,int> bcol2bentry; //maps block column index to dense block entries
383 std::map<int,int>::iterator iter;
384 //Fill the block matrix. We must do this in local index space.
385 //TODO: Optimization: We assume the blocks are fully populated in the point matrix. This means
386 //TODO: on the first point row in the block row, we know that we're hitting new block col indices.
387 //TODO: on other rows, we know the block col indices have all been seen before
388 //int offset;
389 //if (pointMatrix.getIndexBase()) offset = 0;
390 //else offset = 1;
391 for (size_t i=0; i<pointMatrix.getLocalNumRows()/blockSize; i++) {
392 int blkCnt=0; //how many unique block entries encountered so far in current block row
393 for (int j=0; j<blockSize; ++j) {
394 LO rowLid = i*blockSize+j;
395 pointMatrix.getLocalRowView(rowLid,pointColInds,pointVals);
396 for (size_t k=0; k<pointColInds.size(); ++k) {
397 //convert point column to block col
400 if (iter == bcol2bentry.end()) {
401 //new block column
403 blocks[blkCnt].push_back(pointVals[k]);
404 blkCnt++;
405 } else {
406 //block column found previously
407 int littleBlock = iter->second;
408 blocks[littleBlock].push_back(pointVals[k]);
409 }
410 }
411 }
412 // TODO This inserts the blocks one block entry at a time. It is probably more efficient to
413 // TODO store all the blocks in a block row contiguously so they can be inserted with a single call.
414 for (iter=bcol2bentry.begin(); iter != bcol2bentry.end(); ++iter) {
415 LO localBlockCol = iter->first;
416 Scalar *vals = (blocks[iter->second]).getRawPtr();
417 if (std::is_same<typename block_crs_matrix_type::little_block_type::array_layout,Kokkos::LayoutLeft>::value) {
419 for (LO ii=0;ii<blockSize;++ii)
420 for (LO jj=0;jj<blockSize;++jj)
422 Scalar *tmp_vals = tmpBlock.getRawPtr();
423 blockMatrix->replaceLocalValues(i, &localBlockCol, tmp_vals, 1);
424 } else {
426 blockMatrix->replaceLocalValues(i, &localBlockCol, vals, 1);
427 }
428 }
429
430 //Done with block row. Zero everything out.
431 for (int j=0; j<maxBlockEntries; ++j)
432 blocks[j].clear();
433 blkCnt = 0;
434 bcol2bentry.clear();
435 }
436
437 tmpBlock.clear();
438
439 return blockMatrix;
440
441 }
442
443 template<class LO, class GO, class Node>
444 Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
446 {
447 typedef Teuchos::OrdinalTraits<Tpetra::global_size_t> TOT;
448 typedef Tpetra::Map<LO,GO,Node> map_type;
449
450 //calculate mesh GIDs
451 Teuchos::ArrayView<const GO> pointGids = pointMap.getLocalElementList();
452 Teuchos::Array<GO> meshGids;
453 GO indexBase = pointMap.getIndexBase();
454
455 // Use hash table to track whether we've encountered this GID previously. This will happen
456 // when striding through the point DOFs in a block. It should not happen otherwise.
457 // I don't use sort/make unique because I don't want to change the ordering.
458 meshGids.reserve(pointGids.size());
460 for (int i=0; i<pointGids.size(); ++i) {
462 if (hashTable.get(meshGid) == -1) {
463 hashTable.add(meshGid,1); //(key,value)
464 meshGids.push_back(meshGid);
465 }
466 }
467
468 Teuchos::RCP<const map_type> meshMap = Teuchos::rcp( new map_type(TOT::invalid(), meshGids(), 0, pointMap.getComm()) );
469 return meshMap;
470
471 }
472
473
474 template<class LO, class GO, class Node>
475 Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
477 {
478 typedef Teuchos::OrdinalTraits<Tpetra::global_size_t> TOT;
479 typedef Tpetra::Map<LO,GO,Node> map_type;
480
481 //calculate mesh GIDs
482 Teuchos::ArrayView<const GO> blockGids = blockMap.getLocalElementList();
483 Teuchos::Array<GO> pointGids(blockGids.size() * blockSize);
484 GO indexBase = blockMap.getIndexBase();
485
486 for(LO i=0, ct=0; i<(LO)blockGids.size(); i++) {
488 for(LO j=0; j<blockSize; j++, ct++)
490 }
491
492 Teuchos::RCP<const map_type> pointMap = Teuchos::rcp( new map_type(TOT::invalid(), pointGids(), 0, blockMap.getComm()) );
493 return pointMap;
494
495 }
496
497
498 template<class Scalar, class LO, class GO, class Node>
499 Teuchos::RCP<CrsMatrix<Scalar, LO, GO, Node> >
501
502 using Teuchos::Array;
503 using Teuchos::ArrayView;
504 using Teuchos::RCP;
505
506 typedef Tpetra::BlockCrsMatrix<Scalar,LO,GO,Node> block_crs_matrix_type;
507 typedef Tpetra::Map<LO,GO,Node> map_type;
508 typedef Tpetra::CrsMatrix<Scalar, LO,GO,Node> crs_matrix_type;
509
510 using crs_graph_type = typename block_crs_matrix_type::crs_graph_type;
511 using local_graph_device_type = typename crs_matrix_type::local_graph_device_type;
512 using local_matrix_device_type = typename crs_matrix_type::local_matrix_device_type;
513 using row_map_type = typename local_graph_device_type::row_map_type::non_const_type;
514 using entries_type = typename local_graph_device_type::entries_type::non_const_type;
515 using values_type = typename local_matrix_device_type::values_type::non_const_type;
516
517 using row_map_type_const = typename local_graph_device_type::row_map_type;
518 using entries_type_const = typename local_graph_device_type::entries_type;
519
520 using little_block_type = typename block_crs_matrix_type::const_little_block_type;
521 using offset_type = typename row_map_type::non_const_value_type;
522 using column_type = typename entries_type::non_const_value_type;
523
524 using execution_space = typename Node::execution_space;
525 using range_type = Kokkos::RangePolicy<execution_space, LO>;
526
527
528 LO blocksize = blockMatrix.getBlockSize();
529 const offset_type bs2 = blocksize * blocksize;
530 size_t block_nnz = blockMatrix.getLocalNumEntries();
531 size_t point_nnz = block_nnz * bs2;
532
533 // We can get these from the blockMatrix directly
536
537 // We need to generate the row/col Map ourselves.
540
543
544 // Get the last few things
545
546 const crs_graph_type & blockGraph = blockMatrix.getCrsGraph();
547 LO point_rows = (LO) pointRowMap->getLocalNumElements();
548 LO block_rows = (LO) blockRowMap->getLocalNumElements();
549 auto blockValues = blockMatrix.getValuesDevice();
550 auto blockLocalGraph = blockGraph.getLocalGraphDevice();
553
554 // Generate the point matrix rowptr / colind / values
555 row_map_type rowptr("row_map", point_rows+1);
556 entries_type colind("entries", point_nnz);
557 values_type values("values", point_nnz);
558
559 // Pre-fill the rowptr
560 Kokkos::parallel_for("fillRowPtr",range_type(0,block_rows),KOKKOS_LAMBDA(const LO i) {
561 offset_type base = blockRowptr[i];
562 offset_type diff = blockRowptr[i+1] - base;
563 for(LO j=0; j<blocksize; j++) {
564 rowptr[i*blocksize +j] = base*bs2 + j*diff*blocksize;
565 }
566
567 // Fill the last guy, if we're on the final entry
568 if(i==block_rows-1) {
569 rowptr[block_rows*blocksize] = blockRowptr[i+1] * bs2;
570 }
571 });
572
573
574 Kokkos::parallel_for("copyEntriesAndValues",range_type(0,block_rows),KOKKOS_LAMBDA(const LO i) {
575 const offset_type blkBeg = blockRowptr[i];
576 const offset_type numBlocks = blockRowptr[i+1] - blkBeg;
577
578 // For each block in the row...
579 for (offset_type block=0; block < numBlocks; block++) {
581 little_block_type my_block(blockValues.data () + (blkBeg+block) * bs2, blocksize, blocksize);
582
583 // For each entry in the block...
584 for(LO little_row=0; little_row<blocksize; little_row++) {
585 offset_type point_row_offset = rowptr[i*blocksize + little_row];
586 for(LO little_col=0; little_col<blocksize; little_col++) {
589 }
590 }
591
592 }
593 });
594
595 // Generate the matrix
597 pointCrsMatrix->setAllValues(rowptr,colind,values);
598
599 // FUTURE OPTIMIZATION: Directly compute import/export, rather than letting ESFC do it
600 pointCrsMatrix->expertStaticFillComplete(pointDomainMap,pointRangeMap);
601
602 return pointCrsMatrix;
603 }
604
605
606} // namespace Tpetra
607
608//
609// Explicit instantiation macro for blockCrsMatrixWriter (various
610// overloads), writeMatrixStrip, and convertToBlockCrsMatrix.
611//
612// Must be expanded from within the Tpetra namespace!
613//
614#define TPETRA_BLOCKCRSMATRIX_HELPERS_INSTANT(S,LO,GO,NODE) \
615 template void blockCrsMatrixWriter(BlockCrsMatrix<S,LO,GO,NODE> const &A, std::string const &fileName); \
616 template void blockCrsMatrixWriter(BlockCrsMatrix<S,LO,GO,NODE> const &A, std::string const &fileName, Teuchos::ParameterList const &params); \
617 template void blockCrsMatrixWriter(BlockCrsMatrix<S,LO,GO,NODE> const &A, std::ostream &os); \
618 template void blockCrsMatrixWriter(BlockCrsMatrix<S,LO,GO,NODE> const &A, std::ostream &os, Teuchos::ParameterList const &params); \
619 template void writeMatrixStrip(BlockCrsMatrix<S,LO,GO,NODE> const &A, std::ostream &os, Teuchos::ParameterList const &params); \
620 template Teuchos::RCP<BlockCrsMatrix<S, LO, GO, NODE> > convertToBlockCrsMatrix(const CrsMatrix<S, LO, GO, NODE>& pointMatrix, const LO &blockSize); \
621 template Teuchos::RCP<CrsMatrix<S, LO, GO, NODE> > convertToCrsMatrix(const BlockCrsMatrix<S, LO, GO, NODE>& blockMatrix);
622
623//
624// Explicit instantiation macro for createMeshMap / createPointMap
625//
626#define TPETRA_CREATEMESHMAP_INSTANT(LO,GO,NODE) \
627 template Teuchos::RCP<const Map<LO,GO,NODE> > createMeshMap (const LO& blockSize, const Map<LO,GO,NODE>& pointMap); \
628 template Teuchos::RCP<const Map<LO,GO,NODE> > createPointMap (const LO& blockSize, const Map<LO,GO,NODE>& blockMap);
629
630#endif // TPETRA_BLOCKCRSMATRIX_HELPERS_DEF_HPP
Struct that holds views of the contents of a CrsMatrix.
Namespace Tpetra contains the class and methods constituting the Tpetra library.
void writeMatrixStrip(BlockCrsMatrix< Scalar, LO, GO, Node > const &A, std::ostream &os, Teuchos::ParameterList const &params)
Helper function called by blockCrsMatrixWriter.
Teuchos::RCP< const Tpetra::Map< LO, GO, Node > > createPointMap(LO const &blockSize, const Tpetra::Map< LO, GO, Node > &blockMap)
Helper function to generate a point map from a block map (with a given block size) GIDs associated wi...
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 > &params=Teuchos::null)
Nonmember function to create an empty CrsGraph given a row Map and the max number of entries allowed ...
Teuchos::RCP< CrsMatrix< Scalar, LO, GO, Node > > convertToCrsMatrix(const Tpetra::BlockCrsMatrix< Scalar, LO, GO, Node > &blockMatrix)
Non-member constructor that creates a point CrsMatrix from an existing BlockCrsMatrix.
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.
void blockCrsMatrixWriter(BlockCrsMatrix< Scalar, LO, GO, Node > const &A, std::string const &fileName)
Helper function to write a BlockCrsMatrix. Calls the 3-argument version.
@ INSERT
Insert new values that don't currently exist.
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...