Tpetra parallel linear algebra Version of the Day
Loading...
Searching...
No Matches
TpetraExt_TripleMatrixMultiply_def.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Tpetra: Templated Linear Algebra Services Package
4//
5// Copyright 2008 NTESS and the Tpetra contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#ifndef TPETRA_TRIPLEMATRIXMULTIPLY_DEF_HPP
11#define TPETRA_TRIPLEMATRIXMULTIPLY_DEF_HPP
12
13#include "TpetraExt_MatrixMatrix_def.hpp"
14#include "TpetraExt_MatrixMatrix_ExtraKernels_decl.hpp" //for UnmanagedView
15#include "Teuchos_VerboseObject.hpp"
16#include "Teuchos_Array.hpp"
17#include "Tpetra_Util.hpp"
18#include "Tpetra_ConfigDefs.hpp"
19#include "Tpetra_CrsMatrix.hpp"
20#include "TpetraExt_MMHelpers_def.hpp"
21#include "Tpetra_RowMatrixTransposer.hpp"
22#include "Tpetra_ConfigDefs.hpp"
23#include "Tpetra_Map.hpp"
24#include "Tpetra_Export.hpp"
25#include "Tpetra_Import_Util.hpp"
26#include "Tpetra_Import_Util2.hpp"
27#include <algorithm>
28#include <cmath>
29#include "Teuchos_FancyOStream.hpp"
30// #include "KokkosSparse_spgemm.hpp"
31
37/*********************************************************************************************************/
38// Include the architecture-specific kernel partial specializations here
39// NOTE: This needs to be outside all namespaces
40#include "TpetraExt_MatrixMatrix_OpenMP.hpp"
41#include "TpetraExt_MatrixMatrix_Cuda.hpp"
42#include "TpetraExt_MatrixMatrix_HIP.hpp"
43#include "TpetraExt_MatrixMatrix_SYCL.hpp"
44
45namespace Tpetra {
46
47namespace TripleMatrixMultiply {
48
49//
50// This method forms the matrix-matrix product Ac = op(R) * op(A) * op(P), where
51// op(A) == A if transposeA is false,
52// op(A) == A^T if transposeA is true,
53// and similarly for op(R) and op(P).
54//
55template <class Scalar,
56 class LocalOrdinal,
57 class GlobalOrdinal,
58 class Node>
59void MultiplyRAP(
60 const CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& R,
61 bool transposeR,
62 const CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
63 bool transposeA,
64 const CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& P,
65 bool transposeP,
66 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Ac,
67 bool call_FillComplete_on_result,
68 const std::string& label,
69 const Teuchos::RCP<Teuchos::ParameterList>& params) {
70 using Teuchos::null;
71 using Teuchos::RCP;
72 typedef Scalar SC;
73 typedef LocalOrdinal LO;
74 typedef GlobalOrdinal GO;
75 typedef Node NO;
76 typedef CrsMatrix<SC, LO, GO, NO> crs_matrix_type;
77 typedef Import<LO, GO, NO> import_type;
78 typedef Export<LO, GO, NO> export_type;
79 typedef CrsMatrixStruct<SC, LO, GO, NO> crs_matrix_struct_type;
80 typedef Map<LO, GO, NO> map_type;
81 typedef RowMatrixTransposer<SC, LO, GO, NO> transposer_type;
82
83#ifdef HAVE_TPETRA_MMM_TIMINGS
84 std::string prefix_mmm = std::string("TpetraExt ") + label + std::string(": ");
85 using Teuchos::TimeMonitor;
86 RCP<Teuchos::TimeMonitor> MM = rcp(new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string("RAP All Setup"))));
87#endif
88
89 const std::string prefix = "TpetraExt::TripleMatrixMultiply::MultiplyRAP(): ";
90
91 // TEUCHOS_FUNC_TIME_MONITOR_DIFF("My Matrix Mult", mmm_multiply);
92
93 // The input matrices R, A and P must both be fillComplete.
94 TEUCHOS_TEST_FOR_EXCEPTION(!R.isFillComplete(), std::runtime_error, prefix << "Matrix R is not fill complete.");
95 TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), std::runtime_error, prefix << "Matrix A is not fill complete.");
96 TEUCHOS_TEST_FOR_EXCEPTION(!P.isFillComplete(), std::runtime_error, prefix << "Matrix P is not fill complete.");
97
98 // If transposeA is true, then Rprime will be the transpose of R
99 // (computed explicitly via RowMatrixTransposer). Otherwise, Rprime
100 // will just be a pointer to R.
101 RCP<const crs_matrix_type> Rprime = null;
102 // If transposeA is true, then Aprime will be the transpose of A
103 // (computed explicitly via RowMatrixTransposer). Otherwise, Aprime
104 // will just be a pointer to A.
105 RCP<const crs_matrix_type> Aprime = null;
106 // If transposeB is true, then Pprime will be the transpose of P
107 // (computed explicitly via RowMatrixTransposer). Otherwise, Pprime
108 // will just be a pointer to P.
109 RCP<const crs_matrix_type> Pprime = null;
110
111 // Is this a "clean" matrix?
112 //
113 // mfh 27 Sep 2016: Historically, if Epetra_CrsMatrix was neither
114 // locally nor globally indexed, then it was empty. I don't like
115 // this, because the most straightforward implementation presumes
116 // lazy allocation of indices. However, historical precedent
117 // demands that we keep around this predicate as a way to test
118 // whether the matrix is empty.
119 const bool newFlag = !Ac.getGraph()->isLocallyIndexed() && !Ac.getGraph()->isGloballyIndexed();
120
121 using Teuchos::ParameterList;
122 RCP<ParameterList> transposeParams(new ParameterList);
123 transposeParams->set("sort", false);
124
125 if (transposeR && &R != &P) {
126 transposer_type transposer(rcpFromRef(R));
127 Rprime = transposer.createTranspose(transposeParams);
128 } else {
129 Rprime = rcpFromRef(R);
130 }
131
132 if (transposeA) {
133 transposer_type transposer(rcpFromRef(A));
134 Aprime = transposer.createTranspose(transposeParams);
135 } else {
136 Aprime = rcpFromRef(A);
137 }
138
139 if (transposeP) {
140 transposer_type transposer(rcpFromRef(P));
141 Pprime = transposer.createTranspose(transposeParams);
142 } else {
143 Pprime = rcpFromRef(P);
144 }
145
146 // Check size compatibility
147 global_size_t numRCols = R.getDomainMap()->getGlobalNumElements();
148 global_size_t numACols = A.getDomainMap()->getGlobalNumElements();
149 global_size_t numPCols = P.getDomainMap()->getGlobalNumElements();
150 global_size_t Rleft = transposeR ? numRCols : R.getGlobalNumRows();
151 global_size_t Rright = transposeR ? R.getGlobalNumRows() : numRCols;
152 global_size_t Aleft = transposeA ? numACols : A.getGlobalNumRows();
153 global_size_t Aright = transposeA ? A.getGlobalNumRows() : numACols;
154 global_size_t Pleft = transposeP ? numPCols : P.getGlobalNumRows();
155 global_size_t Pright = transposeP ? P.getGlobalNumRows() : numPCols;
156 TEUCHOS_TEST_FOR_EXCEPTION(Rright != Aleft, std::runtime_error,
157 prefix << "ERROR, inner dimensions of op(R) and op(A) "
158 "must match for matrix-matrix product. op(R) is "
159 << Rleft << "x" << Rright << ", op(A) is " << Aleft << "x" << Aright);
160
161 TEUCHOS_TEST_FOR_EXCEPTION(Aright != Pleft, std::runtime_error,
162 prefix << "ERROR, inner dimensions of op(A) and op(P) "
163 "must match for matrix-matrix product. op(A) is "
164 << Aleft << "x" << Aright << ", op(P) is " << Pleft << "x" << Pright);
165
166 // The result matrix Ac must at least have a row-map that reflects the correct
167 // row-size. Don't check the number of columns because rectangular matrices
168 // which were constructed with only one map can still end up having the
169 // correct capacity and dimensions when filled.
170 TEUCHOS_TEST_FOR_EXCEPTION(Rleft > Ac.getGlobalNumRows(), std::runtime_error,
171 prefix << "ERROR, dimensions of result Ac must "
172 "match dimensions of op(R) * op(A) * op(P). Ac has "
173 << Ac.getGlobalNumRows()
174 << " rows, should have at least " << Rleft << std::endl);
175
176 // It doesn't matter whether Ac is already Filled or not. If it is already
177 // Filled, it must have space allocated for the positions that will be
178 // referenced in forming Ac = op(R)*op(A)*op(P). If it doesn't have enough space,
179 // we'll error out later when trying to store result values.
180
181 // CGB: However, matrix must be in active-fill
182 TEUCHOS_TEST_FOR_EXCEPT(Ac.isFillActive() == false);
183
184 // We're going to need to import remotely-owned sections of P if
185 // more than one processor is performing this run, depending on the scenario.
186 int numProcs = P.getComm()->getSize();
187
188 // Declare a couple of structs that will be used to hold views of the data
189 // of R, A and P, to be used for fast access during the matrix-multiplication.
190 crs_matrix_struct_type Rview;
191 crs_matrix_struct_type Aview;
192 crs_matrix_struct_type Pview;
193
194 RCP<const map_type> targetMap_R = Rprime->getRowMap();
195 RCP<const map_type> targetMap_A = Aprime->getRowMap();
196 RCP<const map_type> targetMap_P = Pprime->getRowMap();
197
198#ifdef HAVE_TPETRA_MMM_TIMINGS
199 MM = Teuchos::null;
200 MM = rcp(new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string("RAP All I&X"))));
201#endif
202
203 // Now import any needed remote rows and populate the Aview struct
204 // NOTE: We assert that an import isn't needed --- since we do the transpose
205 // above to handle that.
206 RCP<const import_type> dummyImporter;
207
208 if (!(transposeR && &R == &P))
209 MMdetails::import_and_extract_views(*Rprime, targetMap_R, Rview, dummyImporter, true, label, params);
210
211 MMdetails::import_and_extract_views(*Aprime, targetMap_A, Aview, dummyImporter, true, label, params);
212
213 // We will also need local access to all rows of P that correspond to the
214 // column-map of op(A).
215 if (numProcs > 1)
216 targetMap_P = Aprime->getColMap();
217
218 // Import any needed remote rows and populate the Pview struct.
219 MMdetails::import_and_extract_views(*Pprime, targetMap_P, Pview, Aprime->getGraph()->getImporter(), false, label, params);
220
221 RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Actemp;
222
223 bool needs_final_export = !Pprime->getGraph()->getImporter().is_null();
224 if (needs_final_export)
225 Actemp = rcp(new Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Pprime->getColMap(), 0));
226 else
227 Actemp = rcp(&Ac, false); // don't allow deallocation
228
229#ifdef HAVE_TPETRA_MMM_TIMINGS
230 MM = Teuchos::null;
231 MM = rcp(new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string("RAP All Multiply"))));
232#endif
233
234 // Call the appropriate method to perform the actual multiplication.
235 if (call_FillComplete_on_result && newFlag) {
236 if (transposeR && &R == &P)
237 MMdetails::mult_PT_A_P_newmatrix(Aview, Pview, *Actemp, label, params);
238 else
239 MMdetails::mult_R_A_P_newmatrix(Rview, Aview, Pview, *Actemp, label, params);
240 } else if (call_FillComplete_on_result) {
241 if (transposeR && &R == &P)
242 MMdetails::mult_PT_A_P_reuse(Aview, Pview, *Actemp, label, params);
243 else
244 MMdetails::mult_R_A_P_reuse(Rview, Aview, Pview, *Actemp, label, params);
245 } else {
246 // mfh 27 Sep 2016: Is this the "slow" case? This
247 // "CrsWrapper_CrsMatrix" thing could perhaps be made to support
248 // thread-parallel inserts, but that may take some effort.
249 // CrsWrapper_CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> crsmat(Ac);
250
251 // MMdetails::mult_A_B(Aview, Bview, crsmat, label,params);
252
253 // #ifdef HAVE_TPETRA_MMM_TIMINGS
254 // MM = rcp(new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string("RAP All FillComplete"))));
255 // #endif
256 // if (call_FillComplete_on_result) {
257 // // We'll call FillComplete on the C matrix before we exit, and give it a
258 // // domain-map and a range-map.
259 // // The domain-map will be the domain-map of B, unless
260 // // op(B)==transpose(B), in which case the range-map of B will be used.
261 // // The range-map will be the range-map of A, unless op(A)==transpose(A),
262 // // in which case the domain-map of A will be used.
263 // if (!C.isFillComplete())
264 // C.fillComplete(Bprime->getDomainMap(), Aprime->getRangeMap());
265 // }
266 // Not implemented
267 if (transposeR && &R == &P)
268 MMdetails::mult_PT_A_P_newmatrix(Aview, Pview, *Actemp, label, params);
269 else
270 MMdetails::mult_R_A_P_newmatrix(Rview, Aview, Pview, *Actemp, label, params);
271 }
272
273 if (needs_final_export) {
274#ifdef HAVE_TPETRA_MMM_TIMINGS
275 MM = Teuchos::null;
276 MM = rcp(new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string("RAP exportAndFillComplete"))));
277#endif
278 Teuchos::ParameterList labelList;
279 labelList.set("Timer Label", label);
280 Teuchos::ParameterList& labelList_subList = labelList.sublist("matrixmatrix: kernel params", false);
281 labelList_subList.set("isMatrixMatrix_TransferAndFillComplete", true,
282 "This parameter should be set to true only for MatrixMatrix operations: the optimization in Epetra that was ported to Tpetra does _not_ take into account the possibility that for any given source PID, a particular GID may not exist on the target PID: i.e. a transfer operation. A fix for this general case is in development.");
283
284 RCP<crs_matrix_type> Acprime = rcpFromRef(Ac);
285 if (!params.is_null()) {
286 labelList.setParameters(*params);
287 }
288
289 export_type exporter = export_type(*Pprime->getGraph()->getImporter());
290 Actemp->exportAndFillComplete(Acprime,
291 exporter,
292 Acprime->getDomainMap(),
293 Acprime->getRangeMap(),
294 rcp(&labelList, false));
295 }
296#ifdef HAVE_TPETRA_MMM_STATISTICS
297 printMultiplicationStatistics(Actemp->getGraph()->getExporter(), label + std::string(" RAP MMM"));
298#endif
299}
300
301} // End namespace TripleMatrixMultiply
302
303namespace MMdetails {
304
305// Kernel method for computing the local portion of Ac = R*A*P
306template <class Scalar,
307 class LocalOrdinal,
308 class GlobalOrdinal,
309 class Node>
310void mult_R_A_P_newmatrix(
311 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Rview,
312 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Aview,
313 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Pview,
314 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Ac,
315 const std::string& label,
316 const Teuchos::RCP<Teuchos::ParameterList>& params) {
317 using Teuchos::Array;
318 using Teuchos::ArrayRCP;
319 using Teuchos::ArrayView;
320 using Teuchos::RCP;
321 using Teuchos::rcp;
322
323 // typedef Scalar SC; // unused
324 typedef LocalOrdinal LO;
325 typedef GlobalOrdinal GO;
326 typedef Node NO;
327
328 typedef Import<LO, GO, NO> import_type;
329 typedef Map<LO, GO, NO> map_type;
330
331 // Kokkos typedefs
332 typedef typename map_type::local_map_type local_map_type;
333 typedef typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_matrix_device_type KCRS;
334 typedef typename KCRS::StaticCrsGraphType graph_t;
335 typedef typename graph_t::row_map_type::non_const_type lno_view_t;
336 typedef typename NO::execution_space execution_space;
337 typedef Kokkos::RangePolicy<execution_space, size_t> range_type;
338 typedef Kokkos::View<LO*, typename lno_view_t::array_layout, typename lno_view_t::device_type> lo_view_t;
339
340#ifdef HAVE_TPETRA_MMM_TIMINGS
341 std::string prefix_mmm = std::string("TpetraExt ") + label + std::string(": ");
342 using Teuchos::TimeMonitor;
343 RCP<TimeMonitor> MM = rcp(new TimeMonitor(*(TimeMonitor::getNewTimer(prefix_mmm + std::string("RAP M5 Cmap")))));
344#endif
345 LO LO_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
346
347 // Build the final importer / column map, hash table lookups for Ac
348 RCP<const import_type> Cimport;
349 RCP<const map_type> Ccolmap;
350 RCP<const import_type> Pimport = Pview.origMatrix->getGraph()->getImporter();
351 RCP<const import_type> Iimport = Pview.importMatrix.is_null() ? Teuchos::null : Pview.importMatrix->getGraph()->getImporter();
352 local_map_type Acolmap_local = Aview.colMap->getLocalMap();
353 local_map_type Prowmap_local = Pview.origMatrix->getRowMap()->getLocalMap();
354 local_map_type Irowmap_local;
355 if (!Pview.importMatrix.is_null()) Irowmap_local = Pview.importMatrix->getRowMap()->getLocalMap();
356 local_map_type Pcolmap_local = Pview.origMatrix->getColMap()->getLocalMap();
357 local_map_type Icolmap_local;
358 if (!Pview.importMatrix.is_null()) Icolmap_local = Pview.importMatrix->getColMap()->getLocalMap();
359
360 // mfh 27 Sep 2016: Pcol2Ccol is a table that maps from local column
361 // indices of B, to local column indices of Ac. (B and Ac have the
362 // same number of columns.) The kernel uses this, instead of
363 // copying the entire input matrix B and converting its column
364 // indices to those of C.
365 lo_view_t Pcol2Ccol(Kokkos::ViewAllocateWithoutInitializing("Pcol2Ccol"), Pview.colMap->getLocalNumElements()), Icol2Ccol;
366
367 if (Pview.importMatrix.is_null()) {
368 // mfh 27 Sep 2016: B has no "remotes," so P and C have the same column Map.
369 Cimport = Pimport;
370 Ccolmap = Pview.colMap;
371 const LO colMapSize = static_cast<LO>(Pview.colMap->getLocalNumElements());
372 // Pcol2Ccol is trivial
373 Kokkos::parallel_for(
374 "Tpetra::mult_R_A_P_newmatrix::Pcol2Ccol_fill",
375 Kokkos::RangePolicy<execution_space, LO>(0, colMapSize),
376 KOKKOS_LAMBDA(const LO i) {
377 Pcol2Ccol(i) = i;
378 });
379 } else {
380 // mfh 27 Sep 2016: P has "remotes," so we need to build the
381 // column Map of C, as well as C's Import object (from its domain
382 // Map to its column Map). C's column Map is the union of the
383 // column Maps of (the local part of) P, and the "remote" part of
384 // P. Ditto for the Import. We have optimized this "setUnion"
385 // operation on Import objects and Maps.
386
387 // Choose the right variant of setUnion
388 if (!Pimport.is_null() && !Iimport.is_null()) {
389 Cimport = Pimport->setUnion(*Iimport);
390 } else if (!Pimport.is_null() && Iimport.is_null()) {
391 Cimport = Pimport->setUnion();
392 } else if (Pimport.is_null() && !Iimport.is_null()) {
393 Cimport = Iimport->setUnion();
394 } else {
395 throw std::runtime_error("TpetraExt::RAP status of matrix importers is nonsensical");
396 }
397 Ccolmap = Cimport->getTargetMap();
398
399 // FIXME (mfh 27 Sep 2016) This error check requires an all-reduce
400 // in general. We should get rid of it in order to reduce
401 // communication costs of sparse matrix-matrix multiply.
402 TEUCHOS_TEST_FOR_EXCEPTION(!Cimport->getSourceMap()->isSameAs(*Pview.origMatrix->getDomainMap()),
403 std::runtime_error, "Tpetra::RAP: Import setUnion messed with the DomainMap in an unfortunate way");
404
405 // NOTE: This is not efficient and should be folded into setUnion
406 //
407 // mfh 27 Sep 2016: What the above comment means, is that the
408 // setUnion operation on Import objects could also compute these
409 // local index - to - local index look-up tables.
410 Kokkos::resize(Icol2Ccol, Pview.importMatrix->getColMap()->getLocalNumElements());
411 local_map_type Ccolmap_local = Ccolmap->getLocalMap();
412 Kokkos::parallel_for(
413 "Tpetra::mult_R_A_P_newmatrix::Pcol2Ccol_getGlobalElement", range_type(0, Pview.origMatrix->getColMap()->getLocalNumElements()), KOKKOS_LAMBDA(const LO i) {
414 Pcol2Ccol(i) = Ccolmap_local.getLocalElement(Pcolmap_local.getGlobalElement(i));
415 });
416 Kokkos::parallel_for(
417 "Tpetra::mult_R_A_P_newmatrix::Icol2Ccol_getGlobalElement", range_type(0, Pview.importMatrix->getColMap()->getLocalNumElements()), KOKKOS_LAMBDA(const LO i) {
418 Icol2Ccol(i) = Ccolmap_local.getLocalElement(Icolmap_local.getGlobalElement(i));
419 });
420 }
421
422 // Replace the column map
423 //
424 // mfh 27 Sep 2016: We do this because C was originally created
425 // without a column Map. Now we have its column Map.
426 Ac.replaceColMap(Ccolmap);
427
428 // mfh 27 Sep 2016: Construct tables that map from local column
429 // indices of A, to local row indices of either B_local (the locally
430 // owned part of B), or B_remote (the "imported" remote part of B).
431 //
432 // For column index Aik in row i of A, if the corresponding row of B
433 // exists in the local part of B ("orig") (which I'll call B_local),
434 // then targetMapToOrigRow[Aik] is the local index of that row of B.
435 // Otherwise, targetMapToOrigRow[Aik] is "invalid" (a flag value).
436 //
437 // For column index Aik in row i of A, if the corresponding row of B
438 // exists in the remote part of B ("Import") (which I'll call
439 // B_remote), then targetMapToImportRow[Aik] is the local index of
440 // that row of B. Otherwise, targetMapToOrigRow[Aik] is "invalid"
441 // (a flag value).
442
443 // Run through all the hash table lookups once and for all
444 lo_view_t targetMapToOrigRow(Kokkos::ViewAllocateWithoutInitializing("targetMapToOrigRow"), Aview.colMap->getLocalNumElements());
445 lo_view_t targetMapToImportRow(Kokkos::ViewAllocateWithoutInitializing("targetMapToImportRow"), Aview.colMap->getLocalNumElements());
446 Kokkos::parallel_for(
447 "Tpetra::mult_R_A_P_newmatrix::construct_tables", range_type(Aview.colMap->getMinLocalIndex(), Aview.colMap->getMaxLocalIndex() + 1), KOKKOS_LAMBDA(const LO i) {
448 GO aidx = Acolmap_local.getGlobalElement(i);
449 LO P_LID = Prowmap_local.getLocalElement(aidx);
450 if (P_LID != LO_INVALID) {
451 targetMapToOrigRow(i) = P_LID;
452 targetMapToImportRow(i) = LO_INVALID;
453 } else {
454 LO I_LID = Irowmap_local.getLocalElement(aidx);
455 targetMapToOrigRow(i) = LO_INVALID;
456 targetMapToImportRow(i) = I_LID;
457 }
458 });
459
460 // Call the actual kernel. We'll rely on partial template specialization to call the correct one ---
461 // Either the straight-up Tpetra code (SerialNode) or the KokkosKernels one (other NGP node types)
462 KernelWrappers3<Scalar, LocalOrdinal, GlobalOrdinal, Node, lo_view_t>::
463 mult_R_A_P_newmatrix_kernel_wrapper(Rview, Aview, Pview,
464 targetMapToOrigRow, targetMapToImportRow, Pcol2Ccol, Icol2Ccol,
465 Ac, Cimport, label, params);
466}
467
468// Kernel method for computing the local portion of Ac = R*A*P (reuse mode)
469template <class Scalar,
470 class LocalOrdinal,
471 class GlobalOrdinal,
472 class Node>
473void mult_R_A_P_reuse(
474 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Rview,
475 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Aview,
476 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Pview,
477 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Ac,
478 const std::string& label,
479 const Teuchos::RCP<Teuchos::ParameterList>& params) {
480 using Teuchos::Array;
481 using Teuchos::ArrayRCP;
482 using Teuchos::ArrayView;
483 using Teuchos::RCP;
484 using Teuchos::rcp;
485
486 // typedef Scalar SC; // unused
487 typedef LocalOrdinal LO;
488 typedef GlobalOrdinal GO;
489 typedef Node NO;
490
491 typedef Import<LO, GO, NO> import_type;
492 typedef Map<LO, GO, NO> map_type;
493
494 // Kokkos typedefs
495 typedef typename map_type::local_map_type local_map_type;
496 typedef typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_matrix_device_type KCRS;
497 typedef typename KCRS::StaticCrsGraphType graph_t;
498 typedef typename graph_t::row_map_type::non_const_type lno_view_t;
499 typedef typename NO::execution_space execution_space;
500 typedef Kokkos::RangePolicy<execution_space, size_t> range_type;
501 typedef Kokkos::View<LO*, typename lno_view_t::array_layout, typename lno_view_t::device_type> lo_view_t;
502
503#ifdef HAVE_TPETRA_MMM_TIMINGS
504 std::string prefix_mmm = std::string("TpetraExt ") + label + std::string(": ");
505 using Teuchos::TimeMonitor;
506 RCP<TimeMonitor> MM = rcp(new TimeMonitor(*(TimeMonitor::getNewTimer(prefix_mmm + std::string("RAP M5 Cmap")))));
507#endif
508 LO LO_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
509
510 // Build the final importer / column map, hash table lookups for Ac
511 RCP<const import_type> Cimport = Ac.getGraph()->getImporter();
512 RCP<const map_type> Ccolmap = Ac.getColMap();
513 RCP<const import_type> Pimport = Pview.origMatrix->getGraph()->getImporter();
514 RCP<const import_type> Iimport = Pview.importMatrix.is_null() ? Teuchos::null : Pview.importMatrix->getGraph()->getImporter();
515 local_map_type Acolmap_local = Aview.colMap->getLocalMap();
516 local_map_type Prowmap_local = Pview.origMatrix->getRowMap()->getLocalMap();
517 local_map_type Irowmap_local;
518 if (!Pview.importMatrix.is_null()) Irowmap_local = Pview.importMatrix->getRowMap()->getLocalMap();
519 local_map_type Pcolmap_local = Pview.origMatrix->getColMap()->getLocalMap();
520 local_map_type Icolmap_local;
521 if (!Pview.importMatrix.is_null()) Icolmap_local = Pview.importMatrix->getColMap()->getLocalMap();
522 local_map_type Ccolmap_local = Ccolmap->getLocalMap();
523
524 // Build the final importer / column map, hash table lookups for C
525 lo_view_t Bcol2Ccol(Kokkos::ViewAllocateWithoutInitializing("Bcol2Ccol"), Pview.colMap->getLocalNumElements()), Icol2Ccol;
526 {
527 // Bcol2Col may not be trivial, as Ccolmap is compressed during fillComplete in newmatrix
528 // So, column map of C may be a strict subset of the column map of B
529 Kokkos::parallel_for(
530 range_type(0, Pview.origMatrix->getColMap()->getLocalNumElements()), KOKKOS_LAMBDA(const LO i) {
531 Bcol2Ccol(i) = Ccolmap_local.getLocalElement(Pcolmap_local.getGlobalElement(i));
532 });
533
534 if (!Pview.importMatrix.is_null()) {
535 TEUCHOS_TEST_FOR_EXCEPTION(!Cimport->getSourceMap()->isSameAs(*Pview.origMatrix->getDomainMap()),
536 std::runtime_error, "Tpetra::MMM: Import setUnion messed with the DomainMap in an unfortunate way");
537
538 Kokkos::resize(Icol2Ccol, Pview.importMatrix->getColMap()->getLocalNumElements());
539 Kokkos::parallel_for(
540 range_type(0, Pview.importMatrix->getColMap()->getLocalNumElements()), KOKKOS_LAMBDA(const LO i) {
541 Icol2Ccol(i) = Ccolmap_local.getLocalElement(Icolmap_local.getGlobalElement(i));
542 });
543 }
544 }
545
546 // Run through all the hash table lookups once and for all
547 lo_view_t targetMapToOrigRow(Kokkos::ViewAllocateWithoutInitializing("targetMapToOrigRow"), Aview.colMap->getLocalNumElements());
548 lo_view_t targetMapToImportRow(Kokkos::ViewAllocateWithoutInitializing("targetMapToImportRow"), Aview.colMap->getLocalNumElements());
549 Kokkos::parallel_for(
550 range_type(Aview.colMap->getMinLocalIndex(), Aview.colMap->getMaxLocalIndex() + 1), KOKKOS_LAMBDA(const LO i) {
551 GO aidx = Acolmap_local.getGlobalElement(i);
552 LO B_LID = Prowmap_local.getLocalElement(aidx);
553 if (B_LID != LO_INVALID) {
554 targetMapToOrigRow(i) = B_LID;
555 targetMapToImportRow(i) = LO_INVALID;
556 } else {
557 LO I_LID = Irowmap_local.getLocalElement(aidx);
558 targetMapToOrigRow(i) = LO_INVALID;
559 targetMapToImportRow(i) = I_LID;
560 }
561 });
562
563 // Call the actual kernel. We'll rely on partial template specialization to call the correct one ---
564 // Either the straight-up Tpetra code (SerialNode) or the KokkosKernels one (other NGP node types)
565 KernelWrappers3<Scalar, LocalOrdinal, GlobalOrdinal, Node, lo_view_t>::
566 mult_R_A_P_reuse_kernel_wrapper(Rview, Aview, Pview,
567 targetMapToOrigRow, targetMapToImportRow, Bcol2Ccol, Icol2Ccol,
568 Ac, Cimport, label, params);
569}
570
571// Kernel method for computing the local portion of Ac = R*A*P
572template <class Scalar,
573 class LocalOrdinal,
574 class GlobalOrdinal,
575 class Node>
576void mult_PT_A_P_newmatrix(
577 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Aview,
578 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Pview,
579 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Ac,
580 const std::string& label,
581 const Teuchos::RCP<Teuchos::ParameterList>& params) {
582 using Teuchos::Array;
583 using Teuchos::ArrayRCP;
584 using Teuchos::ArrayView;
585 using Teuchos::RCP;
586 using Teuchos::rcp;
587
588 // typedef Scalar SC; // unused
589 typedef LocalOrdinal LO;
590 typedef GlobalOrdinal GO;
591 typedef Node NO;
592
593 typedef Import<LO, GO, NO> import_type;
594 typedef Map<LO, GO, NO> map_type;
595
596 // Kokkos typedefs
597 typedef typename map_type::local_map_type local_map_type;
598 typedef typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_matrix_device_type KCRS;
599 typedef typename KCRS::StaticCrsGraphType graph_t;
600 typedef typename graph_t::row_map_type::non_const_type lno_view_t;
601 typedef typename NO::execution_space execution_space;
602 typedef Kokkos::RangePolicy<execution_space, size_t> range_type;
603 typedef Kokkos::View<LO*, typename lno_view_t::array_layout, typename lno_view_t::device_type> lo_view_t;
604
605#ifdef HAVE_TPETRA_MMM_TIMINGS
606 std::string prefix_mmm = std::string("TpetraExt ") + label + std::string(": ");
607 using Teuchos::TimeMonitor;
608 RCP<TimeMonitor> MM = rcp(new TimeMonitor(*(TimeMonitor::getNewTimer(prefix_mmm + std::string("PTAP M5 Cmap")))));
609#endif
610 LO LO_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
611
612 // Build the final importer / column map, hash table lookups for Ac
613 RCP<const import_type> Cimport;
614 RCP<const map_type> Ccolmap;
615 RCP<const import_type> Pimport = Pview.origMatrix->getGraph()->getImporter();
616 RCP<const import_type> Iimport = Pview.importMatrix.is_null() ? Teuchos::null : Pview.importMatrix->getGraph()->getImporter();
617 local_map_type Acolmap_local = Aview.colMap->getLocalMap();
618 local_map_type Prowmap_local = Pview.origMatrix->getRowMap()->getLocalMap();
619 local_map_type Irowmap_local;
620 if (!Pview.importMatrix.is_null()) Irowmap_local = Pview.importMatrix->getRowMap()->getLocalMap();
621 local_map_type Pcolmap_local = Pview.origMatrix->getColMap()->getLocalMap();
622 local_map_type Icolmap_local;
623 if (!Pview.importMatrix.is_null()) Icolmap_local = Pview.importMatrix->getColMap()->getLocalMap();
624
625 // mfh 27 Sep 2016: Pcol2Ccol is a table that maps from local column
626 // indices of B, to local column indices of Ac. (B and Ac have the
627 // same number of columns.) The kernel uses this, instead of
628 // copying the entire input matrix B and converting its column
629 // indices to those of C.
630 lo_view_t Pcol2Ccol(Kokkos::ViewAllocateWithoutInitializing("Pcol2Ccol"), Pview.colMap->getLocalNumElements()), Icol2Ccol;
631
632 if (Pview.importMatrix.is_null()) {
633 // mfh 27 Sep 2016: B has no "remotes," so P and C have the same column Map.
634 Cimport = Pimport;
635 Ccolmap = Pview.colMap;
636 const LO colMapSize = static_cast<LO>(Pview.colMap->getLocalNumElements());
637 // Pcol2Ccol is trivial
638 Kokkos::parallel_for(
639 "Tpetra::mult_R_A_P_newmatrix::Pcol2Ccol_fill",
640 Kokkos::RangePolicy<execution_space, LO>(0, colMapSize),
641 KOKKOS_LAMBDA(const LO i) {
642 Pcol2Ccol(i) = i;
643 });
644 } else {
645 // mfh 27 Sep 2016: P has "remotes," so we need to build the
646 // column Map of C, as well as C's Import object (from its domain
647 // Map to its column Map). C's column Map is the union of the
648 // column Maps of (the local part of) P, and the "remote" part of
649 // P. Ditto for the Import. We have optimized this "setUnion"
650 // operation on Import objects and Maps.
651
652 // Choose the right variant of setUnion
653 if (!Pimport.is_null() && !Iimport.is_null()) {
654 Cimport = Pimport->setUnion(*Iimport);
655 } else if (!Pimport.is_null() && Iimport.is_null()) {
656 Cimport = Pimport->setUnion();
657 } else if (Pimport.is_null() && !Iimport.is_null()) {
658 Cimport = Iimport->setUnion();
659 } else {
660 throw std::runtime_error("TpetraExt::RAP status of matrix importers is nonsensical");
661 }
662 Ccolmap = Cimport->getTargetMap();
663
664 // FIXME (mfh 27 Sep 2016) This error check requires an all-reduce
665 // in general. We should get rid of it in order to reduce
666 // communication costs of sparse matrix-matrix multiply.
667 TEUCHOS_TEST_FOR_EXCEPTION(!Cimport->getSourceMap()->isSameAs(*Pview.origMatrix->getDomainMap()),
668 std::runtime_error, "Tpetra::RAP: Import setUnion messed with the DomainMap in an unfortunate way");
669
670 // NOTE: This is not efficient and should be folded into setUnion
671 //
672 // mfh 27 Sep 2016: What the above comment means, is that the
673 // setUnion operation on Import objects could also compute these
674 // local index - to - local index look-up tables.
675 Kokkos::resize(Icol2Ccol, Pview.importMatrix->getColMap()->getLocalNumElements());
676 local_map_type Ccolmap_local = Ccolmap->getLocalMap();
677 Kokkos::parallel_for(
678 "Tpetra::mult_R_A_P_newmatrix::Pcol2Ccol_getGlobalElement", range_type(0, Pview.origMatrix->getColMap()->getLocalNumElements()), KOKKOS_LAMBDA(const LO i) {
679 Pcol2Ccol(i) = Ccolmap_local.getLocalElement(Pcolmap_local.getGlobalElement(i));
680 });
681 Kokkos::parallel_for(
682 "Tpetra::mult_R_A_P_newmatrix::Icol2Ccol_getGlobalElement", range_type(0, Pview.importMatrix->getColMap()->getLocalNumElements()), KOKKOS_LAMBDA(const LO i) {
683 Icol2Ccol(i) = Ccolmap_local.getLocalElement(Icolmap_local.getGlobalElement(i));
684 });
685 }
686
687 // Replace the column map
688 //
689 // mfh 27 Sep 2016: We do this because C was originally created
690 // without a column Map. Now we have its column Map.
691 Ac.replaceColMap(Ccolmap);
692
693 // mfh 27 Sep 2016: Construct tables that map from local column
694 // indices of A, to local row indices of either B_local (the locally
695 // owned part of B), or B_remote (the "imported" remote part of B).
696 //
697 // For column index Aik in row i of A, if the corresponding row of B
698 // exists in the local part of B ("orig") (which I'll call B_local),
699 // then targetMapToOrigRow[Aik] is the local index of that row of B.
700 // Otherwise, targetMapToOrigRow[Aik] is "invalid" (a flag value).
701 //
702 // For column index Aik in row i of A, if the corresponding row of B
703 // exists in the remote part of B ("Import") (which I'll call
704 // B_remote), then targetMapToImportRow[Aik] is the local index of
705 // that row of B. Otherwise, targetMapToOrigRow[Aik] is "invalid"
706 // (a flag value).
707
708 // Run through all the hash table lookups once and for all
709 lo_view_t targetMapToOrigRow(Kokkos::ViewAllocateWithoutInitializing("targetMapToOrigRow"), Aview.colMap->getLocalNumElements());
710 lo_view_t targetMapToImportRow(Kokkos::ViewAllocateWithoutInitializing("targetMapToImportRow"), Aview.colMap->getLocalNumElements());
711
712 Kokkos::parallel_for(
713 "Tpetra::mult_R_A_P_newmatrix::construct_tables", range_type(Aview.colMap->getMinLocalIndex(), Aview.colMap->getMaxLocalIndex() + 1), KOKKOS_LAMBDA(const LO i) {
714 GO aidx = Acolmap_local.getGlobalElement(i);
715 LO P_LID = Prowmap_local.getLocalElement(aidx);
716 if (P_LID != LO_INVALID) {
717 targetMapToOrigRow(i) = P_LID;
718 targetMapToImportRow(i) = LO_INVALID;
719 } else {
720 LO I_LID = Irowmap_local.getLocalElement(aidx);
721 targetMapToOrigRow(i) = LO_INVALID;
722 targetMapToImportRow(i) = I_LID;
723 }
724 });
725
726 // Call the actual kernel. We'll rely on partial template specialization to call the correct one ---
727 // Either the straight-up Tpetra code (SerialNode) or the KokkosKernels one (other NGP node types)
728 KernelWrappers3<Scalar, LocalOrdinal, GlobalOrdinal, Node, lo_view_t>::
729 mult_PT_A_P_newmatrix_kernel_wrapper(Aview, Pview,
730 targetMapToOrigRow, targetMapToImportRow, Pcol2Ccol, Icol2Ccol,
731 Ac, Cimport, label, params);
732}
733
734// Kernel method for computing the local portion of Ac = R*A*P (reuse mode)
735template <class Scalar,
736 class LocalOrdinal,
737 class GlobalOrdinal,
738 class Node>
739void mult_PT_A_P_reuse(
740 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Aview,
741 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Pview,
742 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Ac,
743 const std::string& label,
744 const Teuchos::RCP<Teuchos::ParameterList>& params) {
745 using Teuchos::Array;
746 using Teuchos::ArrayRCP;
747 using Teuchos::ArrayView;
748 using Teuchos::RCP;
749 using Teuchos::rcp;
750
751 // typedef Scalar SC; // unused
752 typedef LocalOrdinal LO;
753 typedef GlobalOrdinal GO;
754 typedef Node NO;
755
756 typedef Import<LO, GO, NO> import_type;
757 typedef Map<LO, GO, NO> map_type;
758
759 // Kokkos typedefs
760 typedef typename map_type::local_map_type local_map_type;
761 typedef typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_matrix_device_type KCRS;
762 typedef typename KCRS::StaticCrsGraphType graph_t;
763 typedef typename graph_t::row_map_type::non_const_type lno_view_t;
764 typedef typename NO::execution_space execution_space;
765 typedef Kokkos::RangePolicy<execution_space, size_t> range_type;
766 typedef Kokkos::View<LO*, typename lno_view_t::array_layout, typename lno_view_t::device_type> lo_view_t;
767
768#ifdef HAVE_TPETRA_MMM_TIMINGS
769 std::string prefix_mmm = std::string("TpetraExt ") + label + std::string(": ");
770 using Teuchos::TimeMonitor;
771 RCP<TimeMonitor> MM = rcp(new TimeMonitor(*(TimeMonitor::getNewTimer(prefix_mmm + std::string("RAP M5 Cmap")))));
772#endif
773 LO LO_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
774
775 // Build the final importer / column map, hash table lookups for Ac
776 RCP<const import_type> Cimport = Ac.getGraph()->getImporter();
777 RCP<const map_type> Ccolmap = Ac.getColMap();
778 RCP<const import_type> Pimport = Pview.origMatrix->getGraph()->getImporter();
779 RCP<const import_type> Iimport = Pview.importMatrix.is_null() ? Teuchos::null : Pview.importMatrix->getGraph()->getImporter();
780 local_map_type Acolmap_local = Aview.colMap->getLocalMap();
781 local_map_type Prowmap_local = Pview.origMatrix->getRowMap()->getLocalMap();
782 local_map_type Irowmap_local;
783 if (!Pview.importMatrix.is_null()) Irowmap_local = Pview.importMatrix->getRowMap()->getLocalMap();
784 local_map_type Pcolmap_local = Pview.origMatrix->getColMap()->getLocalMap();
785 local_map_type Icolmap_local;
786 if (!Pview.importMatrix.is_null()) Icolmap_local = Pview.importMatrix->getColMap()->getLocalMap();
787 local_map_type Ccolmap_local = Ccolmap->getLocalMap();
788
789 // Build the final importer / column map, hash table lookups for C
790 lo_view_t Bcol2Ccol(Kokkos::ViewAllocateWithoutInitializing("Bcol2Ccol"), Pview.colMap->getLocalNumElements()), Icol2Ccol;
791 {
792 // Bcol2Col may not be trivial, as Ccolmap is compressed during fillComplete in newmatrix
793 // So, column map of C may be a strict subset of the column map of B
794 Kokkos::parallel_for(
795 range_type(0, Pview.origMatrix->getColMap()->getLocalNumElements()), KOKKOS_LAMBDA(const LO i) {
796 Bcol2Ccol(i) = Ccolmap_local.getLocalElement(Pcolmap_local.getGlobalElement(i));
797 });
798
799 if (!Pview.importMatrix.is_null()) {
800 TEUCHOS_TEST_FOR_EXCEPTION(!Cimport->getSourceMap()->isSameAs(*Pview.origMatrix->getDomainMap()),
801 std::runtime_error, "Tpetra::MMM: Import setUnion messed with the DomainMap in an unfortunate way");
802
803 Kokkos::resize(Icol2Ccol, Pview.importMatrix->getColMap()->getLocalNumElements());
804 Kokkos::parallel_for(
805 range_type(0, Pview.importMatrix->getColMap()->getLocalNumElements()), KOKKOS_LAMBDA(const LO i) {
806 Icol2Ccol(i) = Ccolmap_local.getLocalElement(Icolmap_local.getGlobalElement(i));
807 });
808 }
809 }
810
811 // Run through all the hash table lookups once and for all
812 lo_view_t targetMapToOrigRow(Kokkos::ViewAllocateWithoutInitializing("targetMapToOrigRow"), Aview.colMap->getLocalNumElements());
813 lo_view_t targetMapToImportRow(Kokkos::ViewAllocateWithoutInitializing("targetMapToImportRow"), Aview.colMap->getLocalNumElements());
814 Kokkos::parallel_for(
815 range_type(Aview.colMap->getMinLocalIndex(), Aview.colMap->getMaxLocalIndex() + 1), KOKKOS_LAMBDA(const LO i) {
816 GO aidx = Acolmap_local.getGlobalElement(i);
817 LO B_LID = Prowmap_local.getLocalElement(aidx);
818 if (B_LID != LO_INVALID) {
819 targetMapToOrigRow(i) = B_LID;
820 targetMapToImportRow(i) = LO_INVALID;
821 } else {
822 LO I_LID = Irowmap_local.getLocalElement(aidx);
823 targetMapToOrigRow(i) = LO_INVALID;
824 targetMapToImportRow(i) = I_LID;
825 }
826 });
827
828 // Call the actual kernel. We'll rely on partial template specialization to call the correct one ---
829 // Either the straight-up Tpetra code (SerialNode) or the KokkosKernels one (other NGP node types)
830 KernelWrappers3<Scalar, LocalOrdinal, GlobalOrdinal, Node, lo_view_t>::
831 mult_PT_A_P_reuse_kernel_wrapper(Aview, Pview,
832 targetMapToOrigRow, targetMapToImportRow, Bcol2Ccol, Icol2Ccol,
833 Ac, Cimport, label, params);
834}
835
836/*********************************************************************************************************/
837// RAP NewMatrix Kernel wrappers (Default non-threaded version)
838// Computes R * A * P -> Ac using classic Gustavson approach
839template <class Scalar,
840 class LocalOrdinal,
841 class GlobalOrdinal,
842 class Node,
843 class LocalOrdinalViewType>
844void KernelWrappers3<Scalar, LocalOrdinal, GlobalOrdinal, Node, LocalOrdinalViewType>::mult_R_A_P_newmatrix_kernel_wrapper(CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Rview,
845 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Aview,
846 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Pview,
847 const LocalOrdinalViewType& Acol2Prow_dev,
848 const LocalOrdinalViewType& Acol2PIrow_dev,
849 const LocalOrdinalViewType& Pcol2Accol_dev,
850 const LocalOrdinalViewType& PIcol2Accol_dev,
851 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Ac,
852 Teuchos::RCP<const Import<LocalOrdinal, GlobalOrdinal, Node> > Acimport,
853 const std::string& label,
854 const Teuchos::RCP<Teuchos::ParameterList>& params) {
855#ifdef HAVE_TPETRA_MMM_TIMINGS
856 std::string prefix_mmm = std::string("TpetraExt ") + label + std::string(": ");
857 using Teuchos::TimeMonitor;
858 Teuchos::RCP<Teuchos::TimeMonitor> MM = rcp(new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string("RAP Newmatrix SerialCore"))));
859#endif
860
861 using Teuchos::Array;
862 using Teuchos::ArrayRCP;
863 using Teuchos::ArrayView;
864 using Teuchos::RCP;
865 using Teuchos::rcp;
866
867 // Lots and lots of typedefs
868 typedef typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_matrix_device_type KCRS;
869 typedef typename KCRS::StaticCrsGraphType graph_t;
870 typedef typename graph_t::row_map_type::const_type c_lno_view_t;
871 typedef typename graph_t::row_map_type::non_const_type lno_view_t;
872 typedef typename graph_t::entries_type::non_const_type lno_nnz_view_t;
873 typedef typename KCRS::values_type::non_const_type scalar_view_t;
874
875 typedef Scalar SC;
876 typedef LocalOrdinal LO;
877 typedef GlobalOrdinal GO;
878 typedef Node NO;
879 typedef Map<LO, GO, NO> map_type;
880 const size_t ST_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
881 const LO LO_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
882 const SC SC_ZERO = Teuchos::ScalarTraits<Scalar>::zero();
883
884 // Sizes
885 RCP<const map_type> Accolmap = Ac.getColMap();
886 size_t m = Rview.origMatrix->getLocalNumRows();
887 size_t n = Accolmap->getLocalNumElements();
888 size_t p_max_nnz_per_row = Pview.origMatrix->getLocalMaxNumRowEntries();
889
890 // Routine runs on host; have to put arguments on host, too
891 auto Acol2Prow = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),
892 Acol2Prow_dev);
893 auto Acol2PIrow = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),
894 Acol2PIrow_dev);
895 auto Pcol2Accol = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),
896 Pcol2Accol_dev);
897 auto PIcol2Accol = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),
898 PIcol2Accol_dev);
899
900 // Grab the Kokkos::SparseCrsMatrices & inner stuff
901 const auto Amat = Aview.origMatrix->getLocalMatrixHost();
902 const auto Pmat = Pview.origMatrix->getLocalMatrixHost();
903 const auto Rmat = Rview.origMatrix->getLocalMatrixHost();
904
905 auto Arowptr = Amat.graph.row_map;
906 auto Prowptr = Pmat.graph.row_map;
907 auto Rrowptr = Rmat.graph.row_map;
908 const auto Acolind = Amat.graph.entries;
909 const auto Pcolind = Pmat.graph.entries;
910 const auto Rcolind = Rmat.graph.entries;
911 const auto Avals = Amat.values;
912 const auto Pvals = Pmat.values;
913 const auto Rvals = Rmat.values;
914
915 typename c_lno_view_t::host_mirror_type::const_type Irowptr;
916 typename lno_nnz_view_t::host_mirror_type Icolind;
917 typename scalar_view_t::host_mirror_type Ivals;
918 if (!Pview.importMatrix.is_null()) {
919 auto lclP = Pview.importMatrix->getLocalMatrixHost();
920 Irowptr = lclP.graph.row_map;
921 Icolind = lclP.graph.entries;
922 Ivals = lclP.values;
923 p_max_nnz_per_row = std::max(p_max_nnz_per_row, Pview.importMatrix->getLocalMaxNumRowEntries());
924 }
925
926#ifdef HAVE_TPETRA_MMM_TIMINGS
927 RCP<TimeMonitor> MM2 = rcp(new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string("RAP Newmatrix SerialCore - Compare"))));
928#endif
929
930 // Classic csr assembly (low memory edition)
931 //
932 // mfh 27 Sep 2016: Ac_estimate_nnz does not promise an upper bound.
933 // The method loops over rows of R, and may resize after processing
934 // each row. Chris Siefert says that this reflects experience in
935 // ML; for the non-threaded case, ML found it faster to spend less
936 // effort on estimation and risk an occasional reallocation.
937 size_t CSR_alloc = std::max(C_estimate_nnz(*Aview.origMatrix, *Pview.origMatrix), n);
938 typename lno_view_t::host_mirror_type Crowptr(Kokkos::ViewAllocateWithoutInitializing("Crowptr"), m + 1);
939 typename lno_nnz_view_t::host_mirror_type Ccolind(Kokkos::ViewAllocateWithoutInitializing("Ccolind"), CSR_alloc);
940 typename scalar_view_t::host_mirror_type Cvals(Kokkos::ViewAllocateWithoutInitializing("Cvals"), CSR_alloc);
941
942 // mfh 27 Sep 2016: The ac_status array is an implementation detail
943 // of the local sparse matrix-matrix multiply routine.
944
945 // The status array will contain the index into colind where this entry was last deposited.
946 // ac_status[i] < nnz - not in the row yet
947 // ac_status[i] >= nnz - this is the entry where you can find the data
948 // We start with this filled with INVALID's indicating that there are no entries yet.
949 // Sadly, this complicates the code due to the fact that size_t's are unsigned.
950 const size_t INVALID = Teuchos::OrdinalTraits<size_t>::invalid();
951 Array<size_t> ac_status(n, ST_INVALID);
952
953 // mfh 27 Sep 2016: Here is the local sparse matrix-matrix multiply
954 // routine. The routine computes Ac := R * A * (P_local + P_remote).
955 //
956 // For column index Aik in row i of A, Acol2Prow[Aik] tells
957 // you whether the corresponding row of P belongs to P_local
958 // ("orig") or P_remote ("Import").
959
960 // For each row of R
961 size_t nnz = 0, nnz_old = 0;
962 for (size_t i = 0; i < m; i++) {
963 // mfh 27 Sep 2016: m is the number of rows in the input matrix R
964 // on the calling process.
965 Crowptr[i] = nnz;
966
967 // mfh 27 Sep 2016: For each entry of R in the current row of R
968 for (size_t kk = Rrowptr[i]; kk < Rrowptr[i + 1]; kk++) {
969 LO k = Rcolind[kk]; // local column index of current entry of R
970 const SC Rik = Rvals[kk]; // value of current entry of R
971 if (Rik == SC_ZERO)
972 continue; // skip explicitly stored zero values in R
973 // For each entry of A in the current row of A
974 for (size_t ll = Arowptr[k]; ll < Arowptr[k + 1]; ll++) {
975 LO l = Acolind[ll]; // local column index of current entry of A
976 const SC Akl = Avals[ll]; // value of current entry of A
977 if (Akl == SC_ZERO)
978 continue; // skip explicitly stored zero values in A
979
980 if (Acol2Prow[l] != LO_INVALID) {
981 // mfh 27 Sep 2016: If the entry of Acol2Prow
982 // corresponding to the current entry of A is populated, then
983 // the corresponding row of P is in P_local (i.e., it lives on
984 // the calling process).
985
986 // Local matrix
987 size_t Pl = Teuchos::as<size_t>(Acol2Prow[l]);
988
989 // mfh 27 Sep 2016: Go through all entries in that row of P_local.
990 for (size_t jj = Prowptr[Pl]; jj < Prowptr[Pl + 1]; jj++) {
991 LO j = Pcolind[jj];
992 LO Acj = Pcol2Accol[j];
993 SC Plj = Pvals[jj];
994
995 if (ac_status[Acj] == INVALID || ac_status[Acj] < nnz_old) {
996#ifdef HAVE_TPETRA_DEBUG
997 // Ac_estimate_nnz() is probably not perfect yet. If this happens, we need to allocate more memory..
998 TEUCHOS_TEST_FOR_EXCEPTION(nnz >= Teuchos::as<size_t>(Ccolind.size()),
999 std::runtime_error,
1000 label << " ERROR, not enough memory allocated for matrix product. Allocated: " << Ccolind.extent(0) << std::endl);
1001#endif
1002 // New entry
1003 ac_status[Acj] = nnz;
1004 Ccolind[nnz] = Acj;
1005 Cvals[nnz] = Rik * Akl * Plj;
1006 nnz++;
1007 } else {
1008 Cvals[ac_status[Acj]] += Rik * Akl * Plj;
1009 }
1010 }
1011 } else {
1012 // mfh 27 Sep 2016: If the entry of Acol2PRow
1013 // corresponding to the current entry of A is NOT populated (has
1014 // a flag "invalid" value), then the corresponding row of P is
1015 // in P_remote (i.e., it does not live on the calling process).
1016
1017 // Remote matrix
1018 size_t Il = Teuchos::as<size_t>(Acol2PIrow[l]);
1019 for (size_t jj = Irowptr[Il]; jj < Irowptr[Il + 1]; jj++) {
1020 LO j = Icolind[jj];
1021 LO Acj = PIcol2Accol[j];
1022 SC Plj = Ivals[jj];
1023
1024 if (ac_status[Acj] == INVALID || ac_status[Acj] < nnz_old) {
1025#ifdef HAVE_TPETRA_DEBUG
1026 // Ac_estimate_nnz() is probably not perfect yet. If this happens, we need to allocate more memory..
1027 TEUCHOS_TEST_FOR_EXCEPTION(nnz >= Teuchos::as<size_t>(Ccolind.size()),
1028 std::runtime_error,
1029 label << " ERROR, not enough memory allocated for matrix product. Allocated: " << Ccolind.extent(0) << std::endl);
1030#endif
1031 // New entry
1032 ac_status[Acj] = nnz;
1033 Ccolind[nnz] = Acj;
1034 Cvals[nnz] = Rik * Akl * Plj;
1035 nnz++;
1036 } else {
1037 Cvals[ac_status[Acj]] += Rik * Akl * Plj;
1038 }
1039 }
1040 }
1041 }
1042 }
1043 // Resize for next pass if needed
1044 if (nnz + n > CSR_alloc) {
1045 CSR_alloc *= 2;
1046 Kokkos::resize(Ccolind, CSR_alloc);
1047 Kokkos::resize(Cvals, CSR_alloc);
1048 }
1049 nnz_old = nnz;
1050 }
1051
1052 Crowptr[m] = nnz;
1053
1054 // Downward resize
1055 Kokkos::resize(Ccolind, nnz);
1056 Kokkos::resize(Cvals, nnz);
1057
1058#ifdef HAVE_TPETRA_MMM_TIMINGS
1059 MM = Teuchos::null;
1060 MM = rcp(new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string("RAP Newmatrix Final Sort"))));
1061#endif
1062 auto Crowptr_dev = Kokkos::create_mirror_view_and_copy(
1063 typename KCRS::device_type(), Crowptr);
1064 auto Ccolind_dev = Kokkos::create_mirror_view_and_copy(
1065 typename KCRS::device_type(), Ccolind);
1066 auto Cvals_dev = Kokkos::create_mirror_view_and_copy(
1067 typename KCRS::device_type(), Cvals);
1068
1069 // Final sort & set of CRS arrays
1070 if (params.is_null() || params->get("sort entries", true))
1071 Import_Util::sortCrsEntries(Crowptr_dev, Ccolind_dev, Cvals_dev);
1072 Ac.setAllValues(Crowptr_dev, Ccolind_dev, Cvals_dev);
1073
1074#ifdef HAVE_TPETRA_MMM_TIMINGS
1075 MM = Teuchos::null;
1076 MM = rcp(new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string("RAP Newmatrix ESFC"))));
1077#endif
1078
1079 // Final FillComplete
1080 //
1081 // mfh 27 Sep 2016: So-called "expert static fill complete" bypasses
1082 // Import (from domain Map to column Map) construction (which costs
1083 // lots of communication) by taking the previously constructed
1084 // Import object. We should be able to do this without interfering
1085 // with the implementation of the local part of sparse matrix-matrix
1086 // multply above.
1087 RCP<Teuchos::ParameterList> labelList = rcp(new Teuchos::ParameterList);
1088 labelList->set("Timer Label", label);
1089 if (!params.is_null()) labelList->set("compute global constants", params->get("compute global constants", true));
1090 RCP<const Export<LO, GO, NO> > dummyExport;
1091 Ac.expertStaticFillComplete(Pview.origMatrix->getDomainMap(),
1092 Rview.origMatrix->getRangeMap(),
1093 Acimport,
1094 dummyExport,
1095 labelList);
1096}
1097
1098/*********************************************************************************************************/
1099// RAP Reuse Kernel wrappers (Default non-threaded version)
1100// Computes R * A * P -> Ac using reuse Gustavson
1101template <class Scalar,
1102 class LocalOrdinal,
1103 class GlobalOrdinal,
1104 class Node,
1105 class LocalOrdinalViewType>
1106void KernelWrappers3<Scalar, LocalOrdinal, GlobalOrdinal, Node, LocalOrdinalViewType>::mult_R_A_P_reuse_kernel_wrapper(CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Rview,
1107 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Aview,
1108 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Pview,
1109 const LocalOrdinalViewType& Acol2Prow_dev,
1110 const LocalOrdinalViewType& Acol2PIrow_dev,
1111 const LocalOrdinalViewType& Pcol2Accol_dev,
1112 const LocalOrdinalViewType& PIcol2Accol_dev,
1113 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Ac,
1114 Teuchos::RCP<const Import<LocalOrdinal, GlobalOrdinal, Node> > Acimport,
1115 const std::string& label,
1116 const Teuchos::RCP<Teuchos::ParameterList>& params) {
1117#ifdef HAVE_TPETRA_MMM_TIMINGS
1118 std::string prefix_mmm = std::string("TpetraExt ") + label + std::string(": ");
1119 using Teuchos::TimeMonitor;
1120 Teuchos::RCP<Teuchos::TimeMonitor> MM = rcp(new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string("RAP Reuse SerialCore"))));
1121#endif
1122
1123 using Teuchos::Array;
1124 using Teuchos::ArrayRCP;
1125 using Teuchos::ArrayView;
1126 using Teuchos::RCP;
1127 using Teuchos::rcp;
1128
1129 // Lots and lots of typedefs
1130 typedef typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_matrix_host_type KCRS;
1131 typedef typename KCRS::StaticCrsGraphType graph_t;
1132 typedef typename graph_t::row_map_type::const_type c_lno_view_t;
1133 typedef typename graph_t::entries_type::non_const_type lno_nnz_view_t;
1134 typedef typename KCRS::values_type::non_const_type scalar_view_t;
1135
1136 typedef Scalar SC;
1137 typedef LocalOrdinal LO;
1138 typedef GlobalOrdinal GO;
1139 typedef Node NO;
1140 typedef Map<LO, GO, NO> map_type;
1141 const size_t ST_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
1142 const LO LO_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
1143 const SC SC_ZERO = Teuchos::ScalarTraits<Scalar>::zero();
1144
1145 // Sizes
1146 RCP<const map_type> Accolmap = Ac.getColMap();
1147 size_t m = Rview.origMatrix->getLocalNumRows();
1148 size_t n = Accolmap->getLocalNumElements();
1149 size_t p_max_nnz_per_row = Pview.origMatrix->getLocalMaxNumRowEntries();
1150
1151 // Routine runs on host; have to put arguments on host, too
1152 auto Acol2Prow = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),
1153 Acol2Prow_dev);
1154 auto Acol2PIrow = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),
1155 Acol2PIrow_dev);
1156 auto Pcol2Accol = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),
1157 Pcol2Accol_dev);
1158 auto PIcol2Accol = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),
1159 PIcol2Accol_dev);
1160
1161 // Grab the Kokkos::SparseCrsMatrices & inner stuff
1162 const KCRS Amat = Aview.origMatrix->getLocalMatrixHost();
1163 const KCRS Pmat = Pview.origMatrix->getLocalMatrixHost();
1164 const KCRS Rmat = Rview.origMatrix->getLocalMatrixHost();
1165 const KCRS Cmat = Ac.getLocalMatrixHost();
1166
1167 c_lno_view_t Arowptr = Amat.graph.row_map, Prowptr = Pmat.graph.row_map, Rrowptr = Rmat.graph.row_map, Crowptr = Cmat.graph.row_map;
1168 const lno_nnz_view_t Acolind = Amat.graph.entries, Pcolind = Pmat.graph.entries, Rcolind = Rmat.graph.entries, Ccolind = Cmat.graph.entries;
1169 const scalar_view_t Avals = Amat.values, Pvals = Pmat.values, Rvals = Rmat.values;
1170 scalar_view_t Cvals = Cmat.values;
1171
1172 c_lno_view_t Irowptr;
1173 lno_nnz_view_t Icolind;
1174 scalar_view_t Ivals;
1175 if (!Pview.importMatrix.is_null()) {
1176 auto lclP = Pview.importMatrix->getLocalMatrixHost();
1177 Irowptr = lclP.graph.row_map;
1178 Icolind = lclP.graph.entries;
1179 Ivals = lclP.values;
1180 p_max_nnz_per_row = std::max(p_max_nnz_per_row, Pview.importMatrix->getLocalMaxNumRowEntries());
1181 }
1182
1183#ifdef HAVE_TPETRA_MMM_TIMINGS
1184 RCP<TimeMonitor> MM2 = rcp(new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string("RAP Reuse SerialCore - Compare"))));
1185#endif
1186
1187 // mfh 27 Sep 2016: The ac_status array is an implementation detail
1188 // of the local sparse matrix-matrix multiply routine.
1189
1190 // The status array will contain the index into colind where this entry was last deposited.
1191 // ac_status[i] < nnz - not in the row yet
1192 // ac_status[i] >= nnz - this is the entry where you can find the data
1193 // We start with this filled with INVALID's indicating that there are no entries yet.
1194 // Sadly, this complicates the code due to the fact that size_t's are unsigned.
1195 Array<size_t> ac_status(n, ST_INVALID);
1196
1197 // mfh 27 Sep 2016: Here is the local sparse matrix-matrix multiply
1198 // routine. The routine computes Ac := R * A * (P_local + P_remote).
1199 //
1200 // For column index Aik in row i of A, Acol2Prow[Aik] tells
1201 // you whether the corresponding row of P belongs to P_local
1202 // ("orig") or P_remote ("Import").
1203
1204 // Necessary until following UVM host accesses are changed - for example Crowptr
1205 // Also probably needed in mult_R_A_P_newmatrix_kernel_wrapper - did not demonstrate this in test failure yet
1206 Kokkos::fence();
1207
1208 // For each row of R
1209 size_t OLD_ip = 0, CSR_ip = 0;
1210 for (size_t i = 0; i < m; i++) {
1211 // First fill the c_status array w/ locations where we're allowed to
1212 // generate nonzeros for this row
1213 OLD_ip = Crowptr[i];
1214 CSR_ip = Crowptr[i + 1];
1215 for (size_t k = OLD_ip; k < CSR_ip; k++) {
1216 ac_status[Ccolind[k]] = k;
1217
1218 // Reset values in the row of C
1219 Cvals[k] = SC_ZERO;
1220 }
1221
1222 // mfh 27 Sep 2016: For each entry of R in the current row of R
1223 for (size_t kk = Rrowptr[i]; kk < Rrowptr[i + 1]; kk++) {
1224 LO k = Rcolind[kk]; // local column index of current entry of R
1225 const SC Rik = Rvals[kk]; // value of current entry of R
1226 if (Rik == SC_ZERO)
1227 continue; // skip explicitly stored zero values in R
1228 // For each entry of A in the current row of A
1229 for (size_t ll = Arowptr[k]; ll < Arowptr[k + 1]; ll++) {
1230 LO l = Acolind[ll]; // local column index of current entry of A
1231 const SC Akl = Avals[ll]; // value of current entry of A
1232 if (Akl == SC_ZERO)
1233 continue; // skip explicitly stored zero values in A
1234
1235 if (Acol2Prow[l] != LO_INVALID) {
1236 // mfh 27 Sep 2016: If the entry of Acol2Prow
1237 // corresponding to the current entry of A is populated, then
1238 // the corresponding row of P is in P_local (i.e., it lives on
1239 // the calling process).
1240
1241 // Local matrix
1242 size_t Pl = Teuchos::as<size_t>(Acol2Prow[l]);
1243
1244 // mfh 27 Sep 2016: Go through all entries in that row of P_local.
1245 for (size_t jj = Prowptr[Pl]; jj < Prowptr[Pl + 1]; jj++) {
1246 LO j = Pcolind[jj];
1247 LO Cij = Pcol2Accol[j];
1248 SC Plj = Pvals[jj];
1249
1250 TEUCHOS_TEST_FOR_EXCEPTION(ac_status[Cij] < OLD_ip || ac_status[Cij] >= CSR_ip,
1251 std::runtime_error, "Trying to insert a new entry (" << i << "," << Cij << ") into a static graph "
1252 << "(c_status = " << ac_status[Cij] << " of [" << OLD_ip << "," << CSR_ip << "))");
1253
1254 Cvals[ac_status[Cij]] += Rik * Akl * Plj;
1255 }
1256 } else {
1257 // mfh 27 Sep 2016: If the entry of Acol2PRow
1258 // corresponding to the current entry of A is NOT populated (has
1259 // a flag "invalid" value), then the corresponding row of P is
1260 // in P_remote (i.e., it does not live on the calling process).
1261
1262 // Remote matrix
1263 size_t Il = Teuchos::as<size_t>(Acol2PIrow[l]);
1264 for (size_t jj = Irowptr[Il]; jj < Irowptr[Il + 1]; jj++) {
1265 LO j = Icolind[jj];
1266 LO Cij = PIcol2Accol[j];
1267 SC Plj = Ivals[jj];
1268
1269 TEUCHOS_TEST_FOR_EXCEPTION(ac_status[Cij] < OLD_ip || ac_status[Cij] >= CSR_ip,
1270 std::runtime_error, "Trying to insert a new entry (" << i << "," << Cij << ") into a static graph "
1271 << "(c_status = " << ac_status[Cij] << " of [" << OLD_ip << "," << CSR_ip << "))");
1272
1273 Cvals[ac_status[Cij]] += Rik * Akl * Plj;
1274 }
1275 }
1276 }
1277 }
1278 }
1279
1280#ifdef HAVE_TPETRA_MMM_TIMINGS
1281 auto MM3 = rcp(new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string("RAP Reuse ESFC"))));
1282#endif
1283
1284 Ac.fillComplete(Ac.getDomainMap(), Ac.getRangeMap());
1285}
1286
1287/*********************************************************************************************************/
1288// PT_A_P NewMatrix Kernel wrappers (Default, general, non-threaded version)
1289// Computes P.T * A * P -> Ac
1290template <class Scalar,
1291 class LocalOrdinal,
1292 class GlobalOrdinal,
1293 class Node,
1294 class LocalOrdinalViewType>
1295void KernelWrappers3<Scalar, LocalOrdinal, GlobalOrdinal, Node, LocalOrdinalViewType>::mult_PT_A_P_newmatrix_kernel_wrapper(CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Aview,
1296 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Pview,
1297 const LocalOrdinalViewType& Acol2Prow,
1298 const LocalOrdinalViewType& Acol2PIrow,
1299 const LocalOrdinalViewType& Pcol2Accol,
1300 const LocalOrdinalViewType& PIcol2Accol,
1301 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Ac,
1302 Teuchos::RCP<const Import<LocalOrdinal, GlobalOrdinal, Node> > Acimport,
1303 const std::string& label,
1304 const Teuchos::RCP<Teuchos::ParameterList>& params) {
1305#ifdef HAVE_TPETRA_MMM_TIMINGS
1306 std::string prefix_mmm = std::string("TpetraExt ") + label + std::string(": ");
1307 using Teuchos::TimeMonitor;
1308 Teuchos::TimeMonitor MM(*TimeMonitor::getNewTimer(prefix_mmm + std::string("PTAP local transpose")));
1309#endif
1310
1311 // We don't need a kernel-level PTAP, we just transpose here
1312 typedef RowMatrixTransposer<Scalar, LocalOrdinal, GlobalOrdinal, Node> transposer_type;
1313 transposer_type transposer(Pview.origMatrix, label + std::string("XP: "));
1314
1315 using Teuchos::ParameterList;
1316 using Teuchos::RCP;
1317 RCP<ParameterList> transposeParams(new ParameterList);
1318 transposeParams->set("sort", false);
1319
1320 if (!params.is_null()) {
1321 transposeParams->set("compute global constants",
1322 params->get("compute global constants: temporaries",
1323 false));
1324 }
1325 RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Ptrans =
1326 transposer.createTransposeLocal(transposeParams);
1327 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node> Rview;
1328 Rview.origMatrix = Ptrans;
1329
1330 mult_R_A_P_newmatrix_kernel_wrapper(Rview, Aview, Pview, Acol2Prow, Acol2PIrow, Pcol2Accol, PIcol2Accol, Ac, Acimport, label, params);
1331}
1332
1333/*********************************************************************************************************/
1334// PT_A_P Reuse Kernel wrappers (Default, general, non-threaded version)
1335// Computes P.T * A * P -> Ac
1336template <class Scalar,
1337 class LocalOrdinal,
1338 class GlobalOrdinal,
1339 class Node,
1340 class LocalOrdinalViewType>
1341void KernelWrappers3<Scalar, LocalOrdinal, GlobalOrdinal, Node, LocalOrdinalViewType>::mult_PT_A_P_reuse_kernel_wrapper(CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Aview,
1342 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Pview,
1343 const LocalOrdinalViewType& Acol2Prow,
1344 const LocalOrdinalViewType& Acol2PIrow,
1345 const LocalOrdinalViewType& Pcol2Accol,
1346 const LocalOrdinalViewType& PIcol2Accol,
1347 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Ac,
1348 Teuchos::RCP<const Import<LocalOrdinal, GlobalOrdinal, Node> > Acimport,
1349 const std::string& label,
1350 const Teuchos::RCP<Teuchos::ParameterList>& params) {
1351#ifdef HAVE_TPETRA_MMM_TIMINGS
1352 std::string prefix_mmm = std::string("TpetraExt ") + label + std::string(": ");
1353 using Teuchos::TimeMonitor;
1354 Teuchos::TimeMonitor MM(*TimeMonitor::getNewTimer(prefix_mmm + std::string("PTAP local transpose")));
1355#endif
1356
1357 // We don't need a kernel-level PTAP, we just transpose here
1358 typedef RowMatrixTransposer<Scalar, LocalOrdinal, GlobalOrdinal, Node> transposer_type;
1359 transposer_type transposer(Pview.origMatrix, label + std::string("XP: "));
1360
1361 using Teuchos::ParameterList;
1362 using Teuchos::RCP;
1363 RCP<ParameterList> transposeParams(new ParameterList);
1364 transposeParams->set("sort", false);
1365
1366 if (!params.is_null()) {
1367 transposeParams->set("compute global constants",
1368 params->get("compute global constants: temporaries",
1369 false));
1370 }
1371 RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Ptrans =
1372 transposer.createTransposeLocal(transposeParams);
1373 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node> Rview;
1374 Rview.origMatrix = Ptrans;
1375
1376 mult_R_A_P_reuse_kernel_wrapper(Rview, Aview, Pview, Acol2Prow, Acol2PIrow, Pcol2Accol, PIcol2Accol, Ac, Acimport, label, params);
1377}
1378
1379/*********************************************************************************************************/
1380// PT_A_P NewMatrix Kernel wrappers (Default non-threaded version)
1381// Computes P.T * A * P -> Ac using a 2-pass algorithm.
1382// This turned out to be slower on SerialNode, but it might still be helpful when going to Kokkos, so I left it in.
1383// Currently, this implementation never gets called.
1384template <class Scalar,
1385 class LocalOrdinal,
1386 class GlobalOrdinal,
1387 class Node>
1388void KernelWrappers3MMM<Scalar, LocalOrdinal, GlobalOrdinal, Node>::mult_PT_A_P_newmatrix_kernel_wrapper_2pass(CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Aview,
1389 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Pview,
1390 const Teuchos::Array<LocalOrdinal>& Acol2PRow,
1391 const Teuchos::Array<LocalOrdinal>& Acol2PRowImport,
1392 const Teuchos::Array<LocalOrdinal>& Pcol2Accol,
1393 const Teuchos::Array<LocalOrdinal>& PIcol2Accol,
1394 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Ac,
1395 Teuchos::RCP<const Import<LocalOrdinal, GlobalOrdinal, Node> > Acimport,
1396 const std::string& label,
1397 const Teuchos::RCP<Teuchos::ParameterList>& params) {
1398#ifdef HAVE_TPETRA_MMM_TIMINGS
1399 std::string prefix_mmm = std::string("TpetraExt ") + label + std::string(": ");
1400 using Teuchos::TimeMonitor;
1401 Teuchos::RCP<Teuchos::TimeMonitor> MM = rcp(new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string("PTAP Newmatrix SerialCore"))));
1402#endif
1403
1404 using Teuchos::Array;
1405 using Teuchos::ArrayRCP;
1406 using Teuchos::ArrayView;
1407 using Teuchos::RCP;
1408 using Teuchos::rcp;
1409
1410 typedef Scalar SC;
1411 typedef LocalOrdinal LO;
1412 typedef GlobalOrdinal GO;
1413 typedef Node NO;
1414 typedef RowMatrixTransposer<SC, LO, GO, NO> transposer_type;
1415 const LO LO_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
1416 const SC SC_ZERO = Teuchos::ScalarTraits<Scalar>::zero();
1417
1418 // number of rows on the process of the fine matrix
1419 // size_t m = Pview.origMatrix->getLocalNumRows();
1420 // number of rows on the process of the coarse matrix
1421 size_t n = Ac.getRowMap()->getLocalNumElements();
1422 LO maxAccol = Ac.getColMap()->getMaxLocalIndex();
1423
1424 // Get Data Pointers
1425 ArrayRCP<size_t> Acrowptr_RCP;
1426 ArrayRCP<LO> Accolind_RCP;
1427 ArrayRCP<SC> Acvals_RCP;
1428
1429 // mfh 27 Sep 2016: get the three CSR arrays
1430 // out of the CrsMatrix. This code computes R * A * (P_local +
1431 // P_remote), where P_local contains the locally owned rows of P,
1432 // and P_remote the (previously Import'ed) remote rows of P.
1433
1434 auto Arowptr = Aview.origMatrix->getLocalRowPtrsHost();
1435 auto Acolind = Aview.origMatrix->getLocalIndicesHost();
1436 auto Avals = Aview.origMatrix->getLocalValuesHost(
1437 Tpetra::Access::ReadOnly);
1438 auto Prowptr = Pview.origMatrix->getLocalRowPtrsHost();
1439 auto Pcolind = Pview.origMatrix->getLocalIndicesHost();
1440 auto Pvals = Pview.origMatrix->getLocalValuesHost(
1441 Tpetra::Access::ReadOnly);
1442 decltype(Prowptr) Irowptr;
1443 decltype(Pcolind) Icolind;
1444 decltype(Pvals) Ivals;
1445
1446 if (!Pview.importMatrix.is_null()) {
1447 Irowptr = Pview.importMatrix->getLocalRowPtrsHost();
1448 Icolind = Pview.importMatrix->getLocalIndicesHost();
1449 Ivals = Pview.importMatrix->getLocalValuesHost(
1450 Tpetra::Access::ReadOnly);
1451 }
1452
1453 // mfh 27 Sep 2016: Remark below "For efficiency" refers to an issue
1454 // where Teuchos::ArrayRCP::operator[] may be slower than
1455 // Teuchos::ArrayView::operator[].
1456
1457 // For efficiency
1458 ArrayView<size_t> Acrowptr;
1459 ArrayView<LO> Accolind;
1460 ArrayView<SC> Acvals;
1461
1463 // In a first pass, determine the graph of Ac.
1465
1467 // Get the graph of Ac. This gets the local transpose of P,
1468 // then loops over R, A, P to get the graph of Ac.
1470
1471#ifdef HAVE_TPETRA_MMM_TIMINGS
1472 MM = rcp(new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string("PTAP local transpose"))));
1473#endif
1474
1476 // Get the local transpose of the graph of P by locally transposing
1477 // all of P
1478
1479 transposer_type transposer(Pview.origMatrix, label + std::string("XP: "));
1480
1481 using Teuchos::ParameterList;
1482 RCP<ParameterList> transposeParams(new ParameterList);
1483 transposeParams->set("sort", false);
1484 if (!params.is_null()) {
1485 transposeParams->set("compute global constants",
1486 params->get("compute global constants: temporaries",
1487 false));
1488 }
1489 RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Ptrans =
1490 transposer.createTransposeLocal(transposeParams);
1491
1492 auto Rrowptr = Ptrans->getLocalRowPtrsHost();
1493 auto Rcolind = Ptrans->getLocalIndicesHost();
1494 auto Rvals = Ptrans->getLocalValuesHost(Tpetra::Access::ReadOnly);
1495
1497 // Construct graph
1498
1499#ifdef HAVE_TPETRA_MMM_TIMINGS
1500 MM = rcp(new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string("PTAP graph"))));
1501#endif
1502
1503 const size_t ST_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
1504 Array<size_t> ac_status(maxAccol + 1, ST_INVALID);
1505
1506 size_t nnz_alloc = std::max(Ac_estimate_nnz(*Aview.origMatrix, *Pview.origMatrix), n);
1507 size_t nnzPerRowA = 100;
1508 if (Aview.origMatrix->getLocalNumEntries() > 0)
1509 nnzPerRowA = Aview.origMatrix->getLocalNumEntries() / Aview.origMatrix->getLocalNumRows();
1510 Acrowptr_RCP.resize(n + 1);
1511 Acrowptr = Acrowptr_RCP();
1512 Accolind_RCP.resize(nnz_alloc);
1513 Accolind = Accolind_RCP();
1514
1515 size_t nnz = 0, nnz_old = 0;
1516 for (size_t i = 0; i < n; i++) {
1517 // mfh 27 Sep 2016: m is the number of rows in the input matrix R
1518 // on the calling process.
1519 Acrowptr[i] = nnz;
1520
1521 // mfh 27 Sep 2016: For each entry of R in the current row of R
1522 for (size_t kk = Rrowptr[i]; kk < Rrowptr[i + 1]; kk++) {
1523 LO k = Rcolind[kk]; // local column index of current entry of R
1524 // For each entry of A in the current row of A
1525 for (size_t ll = Arowptr[k]; ll < Arowptr[k + 1]; ll++) {
1526 LO l = Acolind[ll]; // local column index of current entry of A
1527
1528 if (Acol2PRow[l] != LO_INVALID) {
1529 // mfh 27 Sep 2016: If the entry of Acol2PRow
1530 // corresponding to the current entry of A is populated, then
1531 // the corresponding row of P is in P_local (i.e., it lives on
1532 // the calling process).
1533
1534 // Local matrix
1535 size_t Pl = Teuchos::as<size_t>(Acol2PRow[l]);
1536
1537 // mfh 27 Sep 2016: Go through all entries in that row of P_local.
1538 for (size_t jj = Prowptr[Pl]; jj < Prowptr[Pl + 1]; jj++) {
1539 LO j = Pcolind[jj];
1540 LO Acj = Pcol2Accol[j];
1541
1542 if (ac_status[Acj] == ST_INVALID || ac_status[Acj] < nnz_old) {
1543 // New entry
1544 ac_status[Acj] = nnz;
1545 Accolind[nnz] = Acj;
1546 nnz++;
1547 }
1548 }
1549 } else {
1550 // mfh 27 Sep 2016: If the entry of Acol2PRow
1551 // corresponding to the current entry of A is NOT populated (has
1552 // a flag "invalid" value), then the corresponding row of P is
1553 // in P_remote (i.e., it does not live on the calling process).
1554
1555 // Remote matrix
1556 size_t Il = Teuchos::as<size_t>(Acol2PRowImport[l]);
1557 for (size_t jj = Irowptr[Il]; jj < Irowptr[Il + 1]; jj++) {
1558 LO j = Icolind[jj];
1559 LO Acj = PIcol2Accol[j];
1560
1561 if (ac_status[Acj] == ST_INVALID || ac_status[Acj] < nnz_old) {
1562 // New entry
1563 ac_status[Acj] = nnz;
1564 Accolind[nnz] = Acj;
1565 nnz++;
1566 }
1567 }
1568 }
1569 }
1570 }
1571 // Resize for next pass if needed
1572 // cag: Maybe we can do something more subtle here, and not double
1573 // the size right away.
1574 if (nnz + std::max(5 * nnzPerRowA, n) > nnz_alloc) {
1575 nnz_alloc *= 2;
1576 nnz_alloc = std::max(nnz_alloc, nnz + std::max(5 * nnzPerRowA, n));
1577 Accolind_RCP.resize(nnz_alloc);
1578 Accolind = Accolind_RCP();
1579 Acvals_RCP.resize(nnz_alloc);
1580 Acvals = Acvals_RCP();
1581 }
1582 nnz_old = nnz;
1583 }
1584 Acrowptr[n] = nnz;
1585
1586 // Downward resize
1587 Accolind_RCP.resize(nnz);
1588 Accolind = Accolind_RCP();
1589
1590 // Allocate Acvals
1591 Acvals_RCP.resize(nnz, SC_ZERO);
1592 Acvals = Acvals_RCP();
1593
1595 // In a second pass, enter the values into Acvals.
1597
1598#ifdef HAVE_TPETRA_MMM_TIMINGS
1599 MM = rcp(new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string("PTAP Newmatrix Fill Matrix"))));
1600#endif
1601
1602 for (size_t k = 0; k < n; k++) {
1603 for (size_t ii = Prowptr[k]; ii < Prowptr[k + 1]; ii++) {
1604 LO i = Pcolind[ii];
1605 const SC Pki = Pvals[ii];
1606 for (size_t ll = Arowptr[k]; ll < Arowptr[k + 1]; ll++) {
1607 LO l = Acolind[ll];
1608 const SC Akl = Avals[ll];
1609 if (Akl == 0.)
1610 continue;
1611 if (Acol2PRow[l] != LO_INVALID) {
1612 // mfh 27 Sep 2016: If the entry of Acol2PRow
1613 // corresponding to the current entry of A is populated, then
1614 // the corresponding row of P is in P_local (i.e., it lives on
1615 // the calling process).
1616
1617 // Local matrix
1618 size_t Pl = Teuchos::as<size_t>(Acol2PRow[l]);
1619 for (size_t jj = Prowptr[Pl]; jj < Prowptr[Pl + 1]; jj++) {
1620 LO j = Pcolind[jj];
1621 LO Acj = Pcol2Accol[j];
1622 size_t pp;
1623 for (pp = Acrowptr[i]; pp < Acrowptr[i + 1]; pp++)
1624 if (Accolind[pp] == Acj)
1625 break;
1626 // TEUCHOS_TEST_FOR_EXCEPTION(Accolind[pp] != Acj,
1627 // std::runtime_error, "problem with Ac column indices");
1628 Acvals[pp] += Pki * Akl * Pvals[jj];
1629 }
1630 } else {
1631 // mfh 27 Sep 2016: If the entry of Acol2PRow
1632 // corresponding to the current entry of A NOT populated (has
1633 // a flag "invalid" value), then the corresponding row of P is
1634 // in P_remote (i.e., it does not live on the calling process).
1635
1636 // Remote matrix
1637 size_t Il = Teuchos::as<size_t>(Acol2PRowImport[l]);
1638 for (size_t jj = Irowptr[Il]; jj < Irowptr[Il + 1]; jj++) {
1639 LO j = Icolind[jj];
1640 LO Acj = PIcol2Accol[j];
1641 size_t pp;
1642 for (pp = Acrowptr[i]; pp < Acrowptr[i + 1]; pp++)
1643 if (Accolind[pp] == Acj)
1644 break;
1645 // TEUCHOS_TEST_FOR_EXCEPTION(Accolind[pp] != Acj,
1646 // std::runtime_error, "problem with Ac column indices");
1647 Acvals[pp] += Pki * Akl * Ivals[jj];
1648 }
1649 }
1650 }
1651 }
1652 }
1653
1654#ifdef HAVE_TPETRA_MMM_TIMINGS
1655 MM = rcp(new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string("PTAP sort"))));
1656#endif
1657
1658 // Final sort & set of CRS arrays
1659 //
1660 // TODO (mfh 27 Sep 2016) Will the thread-parallel "local" sparse
1661 // matrix-matrix multiply routine sort the entries for us?
1662 Import_Util::sortCrsEntries(Acrowptr_RCP(), Accolind_RCP(), Acvals_RCP());
1663
1664 // mfh 27 Sep 2016: This just sets pointers.
1665 Ac.setAllValues(Acrowptr_RCP, Accolind_RCP, Acvals_RCP);
1666
1667#ifdef HAVE_TPETRA_MMM_TIMINGS
1668 MM = rcp(new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string("PTAP Newmatrix ESFC"))));
1669#endif
1670
1671 // Final FillComplete
1672 //
1673 // mfh 27 Sep 2016: So-called "expert static fill complete" bypasses
1674 // Import (from domain Map to column Map) construction (which costs
1675 // lots of communication) by taking the previously constructed
1676 // Import object. We should be able to do this without interfering
1677 // with the implementation of the local part of sparse matrix-matrix
1678 // multply above.
1679 RCP<Teuchos::ParameterList> labelList = rcp(new Teuchos::ParameterList);
1680 labelList->set("Timer Label", label);
1681 // labelList->set("Sort column Map ghost GIDs")
1682 if (!params.is_null()) labelList->set("compute global constants", params->get("compute global constants", true));
1683 RCP<const Export<LO, GO, NO> > dummyExport;
1684 Ac.expertStaticFillComplete(Pview.origMatrix->getDomainMap(),
1685 Pview.origMatrix->getDomainMap(),
1686 Acimport,
1687 dummyExport, labelList);
1688}
1689
1690} // namespace MMdetails
1691
1692} // End namespace Tpetra
1693//
1694// Explicit instantiation macro
1695//
1696// Must be expanded from within the Tpetra namespace!
1697//
1698
1699#define TPETRA_TRIPLEMATRIXMULTIPLY_INSTANT(SCALAR, LO, GO, NODE) \
1700 \
1701 template void TripleMatrixMultiply::MultiplyRAP( \
1702 const CrsMatrix<SCALAR, LO, GO, NODE>& R, \
1703 bool transposeR, \
1704 const CrsMatrix<SCALAR, LO, GO, NODE>& A, \
1705 bool transposeA, \
1706 const CrsMatrix<SCALAR, LO, GO, NODE>& P, \
1707 bool transposeP, \
1708 CrsMatrix<SCALAR, LO, GO, NODE>& Ac, \
1709 bool call_FillComplete_on_result, \
1710 const std::string& label, \
1711 const Teuchos::RCP<Teuchos::ParameterList>& params);
1712
1713#endif // TPETRA_TRIPLEMATRIXMULTIPLY_DEF_HPP
Tpetra_Import_Util2.hpp
Utility functions for packing and unpacking sparse matrix entries.
Tpetra_Import_Util.hpp
Internal functions and macros designed for use with Tpetra::Import and Tpetra::Export objects.
Tpetra_Util.hpp
Stand-alone utility functions and macros.
TpetraExt_MMHelpers_def.hpp
TpetraExt_MatrixMatrix_def.hpp
Tpetra::CrsMatrix::local_matrix_device_type
KokkosSparse::CrsMatrix< impl_scalar_type, local_ordinal_type, device_type, void, typename local_graph_device_type::size_type > local_matrix_device_type
The specialization of Kokkos::CrsMatrix that represents the part of the sparse matrix on each MPI pro...
Definition Tpetra_CrsMatrix_decl.hpp:472
Tpetra::CrsMatrixStruct
Struct that holds views of the contents of a CrsMatrix.
Definition TpetraExt_MMHelpers_decl.hpp:36
Tpetra::TripleMatrixMultiply::MultiplyRAP
void MultiplyRAP(const CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &R, bool transposeR, const CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, bool transposeA, const CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &P, bool transposeP, CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Ac, bool call_FillComplete_on_result=true, const std::string &label=std::string(), const Teuchos::RCP< Teuchos::ParameterList > &params=Teuchos::null)
Sparse matrix-matrix multiply.
Definition TpetraExt_TripleMatrixMultiply_def.hpp:59
Tpetra
Namespace Tpetra contains the class and methods constituting the Tpetra library.
Generated by doxygen 1.9.8