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