41 #ifndef TPETRA_TRIPLEMATRIXMULTIPLY_DEF_HPP
42 #define TPETRA_TRIPLEMATRIXMULTIPLY_DEF_HPP
45 #include "TpetraExt_MatrixMatrix_ExtraKernels_decl.hpp"
46 #include "Teuchos_VerboseObject.hpp"
47 #include "Teuchos_Array.hpp"
49 #include "Tpetra_ConfigDefs.hpp"
50 #include "Tpetra_CrsMatrix.hpp"
52 #include "Tpetra_RowMatrixTransposer.hpp"
53 #include "Tpetra_ConfigDefs.hpp"
54 #include "Tpetra_Map.hpp"
55 #include "Tpetra_Export.hpp"
60 #include "Teuchos_FancyOStream.hpp"
74 #include "TpetraExt_MatrixMatrix_OpenMP.hpp"
75 #include "TpetraExt_MatrixMatrix_Cuda.hpp"
79 namespace TripleMatrixMultiply{
87 template <
class Scalar,
99 bool call_FillComplete_on_result,
100 const std::string& label,
101 const Teuchos::RCP<Teuchos::ParameterList>& params)
106 typedef LocalOrdinal LO;
107 typedef GlobalOrdinal GO;
116 #ifdef HAVE_TPETRA_MMM_TIMINGS
117 std::string prefix_mmm = std::string(
"TpetraExt ") + label + std::string(
": ");
118 using Teuchos::TimeMonitor;
119 RCP<Teuchos::TimeMonitor> MM = rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string(
"RAP All Setup"))));
122 const std::string prefix =
"TpetraExt::TripleMatrixMultiply::MultiplyRAP(): ";
127 TEUCHOS_TEST_FOR_EXCEPTION(!R.
isFillComplete(), std::runtime_error, prefix <<
"Matrix R is not fill complete.");
128 TEUCHOS_TEST_FOR_EXCEPTION(!A.
isFillComplete(), std::runtime_error, prefix <<
"Matrix A is not fill complete.");
129 TEUCHOS_TEST_FOR_EXCEPTION(!P.
isFillComplete(), std::runtime_error, prefix <<
"Matrix P is not fill complete.");
134 RCP<const crs_matrix_type> Rprime = null;
138 RCP<const crs_matrix_type> Aprime = null;
142 RCP<const crs_matrix_type> Pprime = null;
152 const bool newFlag = !Ac.
getGraph()->isLocallyIndexed() && !Ac.
getGraph()->isGloballyIndexed();
154 using Teuchos::ParameterList;
155 RCP<ParameterList> transposeParams (
new ParameterList);
156 transposeParams->set (
"sort",
false);
158 if (transposeR && &R != &P) {
159 transposer_type transposer(rcpFromRef (R));
160 Rprime = transposer.createTranspose (transposeParams);
162 Rprime = rcpFromRef(R);
166 transposer_type transposer(rcpFromRef (A));
167 Aprime = transposer.createTranspose (transposeParams);
169 Aprime = rcpFromRef(A);
173 transposer_type transposer(rcpFromRef (P));
174 Pprime = transposer.createTranspose (transposeParams);
176 Pprime = rcpFromRef(P);
189 TEUCHOS_TEST_FOR_EXCEPTION(Rright != Aleft, std::runtime_error,
190 prefix <<
"ERROR, inner dimensions of op(R) and op(A) "
191 "must match for matrix-matrix product. op(R) is "
192 << Rleft <<
"x" << Rright <<
", op(A) is "<< Aleft <<
"x" << Aright);
194 TEUCHOS_TEST_FOR_EXCEPTION(Aright != Pleft, std::runtime_error,
195 prefix <<
"ERROR, inner dimensions of op(A) and op(P) "
196 "must match for matrix-matrix product. op(A) is "
197 << Aleft <<
"x" << Aright <<
", op(P) is "<< Pleft <<
"x" << Pright);
203 TEUCHOS_TEST_FOR_EXCEPTION(Rleft > Ac.
getGlobalNumRows(), std::runtime_error,
204 prefix <<
"ERROR, dimensions of result Ac must "
205 "match dimensions of op(R) * op(A) * op(P). Ac has " << Ac.
getGlobalNumRows()
206 <<
" rows, should have at least " << Rleft << std::endl);
218 int numProcs = P.
getComm()->getSize();
222 crs_matrix_struct_type Rview;
223 crs_matrix_struct_type Aview;
224 crs_matrix_struct_type Pview;
226 RCP<const map_type> targetMap_R = Rprime->getRowMap();
227 RCP<const map_type> targetMap_A = Aprime->getRowMap();
228 RCP<const map_type> targetMap_P = Pprime->getRowMap();
230 #ifdef HAVE_TPETRA_MMM_TIMINGS
231 MM = rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string(
"RAP All I&X"))));
237 RCP<const import_type> dummyImporter;
239 if (!(transposeR && &R == &P))
240 MMdetails::import_and_extract_views(*Rprime, targetMap_R, Rview, dummyImporter,
true, label, params);
242 MMdetails::import_and_extract_views(*Aprime, targetMap_A, Aview, dummyImporter,
true, label, params);
247 targetMap_P = Aprime->getColMap();
250 MMdetails::import_and_extract_views(*Pprime, targetMap_P, Pview, Aprime->getGraph()->getImporter(),
false, label, params);
253 RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Actemp;
255 bool needs_final_export = !Pprime->getGraph()->getImporter().is_null();
256 if (needs_final_export)
259 Actemp = rcp(&Ac,
false);
261 #ifdef HAVE_TPETRA_MMM_TIMINGS
262 MM = rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string(
"RAP All Multiply"))));
266 if (call_FillComplete_on_result && newFlag) {
267 if (transposeR && &R == &P)
268 MMdetails::mult_PT_A_P_newmatrix(Aview, Pview, *Actemp, label, params);
270 MMdetails::mult_R_A_P_newmatrix(Rview, Aview, Pview, *Actemp, label, params);
271 }
else if (call_FillComplete_on_result) {
272 if (transposeR && &R == &P)
273 MMdetails::mult_PT_A_P_reuse(Aview, Pview, *Actemp, label, params);
275 MMdetails::mult_R_A_P_reuse(Rview, Aview, Pview, *Actemp, label, params);
298 if (transposeR && &R == &P)
299 MMdetails::mult_PT_A_P_newmatrix(Aview, Pview, *Actemp, label, params);
301 MMdetails::mult_R_A_P_newmatrix(Rview, Aview, Pview, *Actemp, label, params);
304 if (needs_final_export) {
305 #ifdef HAVE_TPETRA_MMM_TIMINGS
306 MM = rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string(
"RAP exportAndFillComplete"))));
308 Teuchos::ParameterList labelList;
309 labelList.set(
"Timer Label", label);
310 Teuchos::ParameterList& labelList_subList = labelList.sublist(
"matrixmatrix: kernel params",
false);
312 RCP<crs_matrix_type> Acprime = rcpFromRef(Ac);
314 bool overrideAllreduce =
false;
316 if(!params.is_null()) {
317 Teuchos::ParameterList& params_sublist = params->sublist(
"matrixmatrix: kernel params",
false);
319 mm_optimization_core_count = params_sublist.get(
"MM_TAFC_OptimizationCoreCount",mm_optimization_core_count);
320 int mm_optimization_core_count2 = params->get(
"MM_TAFC_OptimizationCoreCount",mm_optimization_core_count);
321 if(mm_optimization_core_count2<mm_optimization_core_count) mm_optimization_core_count=mm_optimization_core_count2;
322 isMM = params_sublist.get(
"isMatrixMatrix_TransferAndFillComplete",
false);
323 overrideAllreduce = params_sublist.get(
"MM_TAFC_OverrideAllreduceCheck",
false);
325 labelList.set(
"compute global constants",params->get(
"compute global constants",
true));
327 labelList_subList.set(
"MM_TAFC_OptimizationCoreCount",mm_optimization_core_count,
"Core Count above which the optimized neighbor discovery is used");
329 labelList_subList.set(
"isMatrixMatrix_TransferAndFillComplete",isMM,
330 "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.");
331 labelList_subList.set(
"MM_TAFC_OverrideAllreduceCheck",overrideAllreduce);
333 export_type exporter = export_type(*Pprime->getGraph()->getImporter());
334 Actemp->exportAndFillComplete(Acprime,
336 Acprime->getDomainMap(),
337 Acprime->getRangeMap(),
338 rcp(&labelList,
false));
341 #ifdef HAVE_TPETRA_MMM_STATISTICS
342 printMultiplicationStatistics(Actemp->getGraph()->getExporter(), label+std::string(
" RAP MMM"));
353 template<
class Scalar,
357 void mult_R_A_P_newmatrix(
358 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Rview,
359 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Aview,
360 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Pview,
361 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Ac,
362 const std::string& label,
363 const Teuchos::RCP<Teuchos::ParameterList>& params)
365 using Teuchos::Array;
366 using Teuchos::ArrayRCP;
367 using Teuchos::ArrayView;
372 typedef LocalOrdinal LO;
373 typedef GlobalOrdinal GO;
376 typedef Import<LO,GO,NO> import_type;
377 typedef Map<LO,GO,NO> map_type;
380 typedef typename map_type::local_map_type local_map_type;
382 typedef typename KCRS::StaticCrsGraphType graph_t;
383 typedef typename graph_t::row_map_type::non_const_type lno_view_t;
384 typedef typename NO::execution_space execution_space;
385 typedef Kokkos::RangePolicy<execution_space, size_t> range_type;
386 typedef Kokkos::View<LO*, typename lno_view_t::array_layout, typename lno_view_t::device_type> lo_view_t;
388 #ifdef HAVE_TPETRA_MMM_TIMINGS
389 std::string prefix_mmm = std::string(
"TpetraExt ") + label + std::string(
": ");
390 using Teuchos::TimeMonitor;
391 RCP<TimeMonitor> MM = rcp(
new TimeMonitor(*(TimeMonitor::getNewTimer(prefix_mmm + std::string(
"RAP M5 Cmap")))));
393 LO LO_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
396 RCP<const import_type> Cimport;
397 RCP<const map_type> Ccolmap;
398 RCP<const import_type> Pimport = Pview.origMatrix->getGraph()->getImporter();
399 RCP<const import_type> Iimport = Pview.importMatrix.is_null() ? Teuchos::null : Pview.importMatrix->getGraph()->getImporter();
400 local_map_type Acolmap_local = Aview.colMap->getLocalMap();
401 local_map_type Prowmap_local = Pview.origMatrix->getRowMap()->getLocalMap();
402 local_map_type Irowmap_local;
if(!Pview.importMatrix.is_null()) Irowmap_local = Pview.importMatrix->getRowMap()->getLocalMap();
403 local_map_type Pcolmap_local = Pview.origMatrix->getColMap()->getLocalMap();
404 local_map_type Icolmap_local;
if(!Pview.importMatrix.is_null()) Icolmap_local = Pview.importMatrix->getColMap()->getLocalMap();
412 lo_view_t Pcol2Ccol(Kokkos::ViewAllocateWithoutInitializing(
"Pcol2Ccol"),Pview.colMap->getNodeNumElements()), Icol2Ccol;
414 if (Pview.importMatrix.is_null()) {
417 Ccolmap = Pview.colMap;
418 const LO colMapSize =
static_cast<LO
>(Pview.colMap->getNodeNumElements());
420 Kokkos::parallel_for(
"Tpetra::mult_R_A_P_newmatrix::Pcol2Ccol_fill",
421 Kokkos::RangePolicy<execution_space, LO>(0, colMapSize),
422 KOKKOS_LAMBDA(
const LO i) {
435 if (!Pimport.is_null() && !Iimport.is_null()) {
436 Cimport = Pimport->setUnion(*Iimport);
438 else if (!Pimport.is_null() && Iimport.is_null()) {
439 Cimport = Pimport->setUnion();
441 else if (Pimport.is_null() && !Iimport.is_null()) {
442 Cimport = Iimport->setUnion();
445 throw std::runtime_error(
"TpetraExt::RAP status of matrix importers is nonsensical");
447 Ccolmap = Cimport->getTargetMap();
452 TEUCHOS_TEST_FOR_EXCEPTION(!Cimport->getSourceMap()->isSameAs(*Pview.origMatrix->getDomainMap()),
453 std::runtime_error,
"Tpetra::RAP: Import setUnion messed with the DomainMap in an unfortunate way");
460 Kokkos::resize(Icol2Ccol,Pview.importMatrix->getColMap()->getNodeNumElements());
461 local_map_type Ccolmap_local = Ccolmap->getLocalMap();
462 Kokkos::parallel_for(
"Tpetra::mult_R_A_P_newmatrix::Pcol2Ccol_getGlobalElement",range_type(0,Pview.origMatrix->getColMap()->getNodeNumElements()),KOKKOS_LAMBDA(
const LO i) {
463 Pcol2Ccol(i) = Ccolmap_local.getLocalElement(Pcolmap_local.getGlobalElement(i));
465 Kokkos::parallel_for(
"Tpetra::mult_R_A_P_newmatrix::Icol2Ccol_getGlobalElement",range_type(0,Pview.importMatrix->getColMap()->getNodeNumElements()),KOKKOS_LAMBDA(
const LO i) {
466 Icol2Ccol(i) = Ccolmap_local.getLocalElement(Icolmap_local.getGlobalElement(i));
475 Ac.replaceColMap(Ccolmap);
493 lo_view_t targetMapToOrigRow(Kokkos::ViewAllocateWithoutInitializing(
"targetMapToOrigRow"),Aview.colMap->getNodeNumElements());
494 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;
503 LO I_LID = Irowmap_local.getLocalElement(aidx);
504 targetMapToOrigRow(i) = LO_INVALID;
505 targetMapToImportRow(i) = I_LID;
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);
519 template<
class Scalar,
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)
531 using Teuchos::Array;
532 using Teuchos::ArrayRCP;
533 using Teuchos::ArrayView;
538 typedef LocalOrdinal LO;
539 typedef GlobalOrdinal GO;
542 typedef Import<LO,GO,NO> import_type;
543 typedef Map<LO,GO,NO> map_type;
546 typedef typename map_type::local_map_type local_map_type;
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;
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")))));
559 LO LO_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
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();
574 lo_view_t Bcol2Ccol(Kokkos::ViewAllocateWithoutInitializing(
"Bcol2Ccol"),Pview.colMap->getNodeNumElements()), Icol2Ccol;
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));
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");
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));
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;
603 LO I_LID = Irowmap_local.getLocalElement(aidx);
604 targetMapToOrigRow(i) = LO_INVALID;
605 targetMapToImportRow(i) = I_LID;
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);
620 template<
class Scalar,
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)
631 using Teuchos::Array;
632 using Teuchos::ArrayRCP;
633 using Teuchos::ArrayView;
638 typedef LocalOrdinal LO;
639 typedef GlobalOrdinal GO;
642 typedef Import<LO,GO,NO> import_type;
643 typedef Map<LO,GO,NO> map_type;
646 typedef typename map_type::local_map_type local_map_type;
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;
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")))));
659 LO LO_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
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();
678 lo_view_t Pcol2Ccol(Kokkos::ViewAllocateWithoutInitializing(
"Pcol2Ccol"),Pview.colMap->getNodeNumElements()), Icol2Ccol;
680 if (Pview.importMatrix.is_null()) {
683 Ccolmap = Pview.colMap;
684 const LO colMapSize =
static_cast<LO
>(Pview.colMap->getNodeNumElements());
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) {
701 if (!Pimport.is_null() && !Iimport.is_null()) {
702 Cimport = Pimport->setUnion(*Iimport);
704 else if (!Pimport.is_null() && Iimport.is_null()) {
705 Cimport = Pimport->setUnion();
707 else if (Pimport.is_null() && !Iimport.is_null()) {
708 Cimport = Iimport->setUnion();
711 throw std::runtime_error(
"TpetraExt::RAP status of matrix importers is nonsensical");
713 Ccolmap = Cimport->getTargetMap();
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");
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));
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));
741 Ac.replaceColMap(Ccolmap);
759 lo_view_t targetMapToOrigRow(Kokkos::ViewAllocateWithoutInitializing(
"targetMapToOrigRow"),Aview.colMap->getNodeNumElements());
760 lo_view_t targetMapToImportRow(Kokkos::ViewAllocateWithoutInitializing(
"targetMapToImportRow"),Aview.colMap->getNodeNumElements());
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;
769 LO I_LID = Irowmap_local.getLocalElement(aidx);
770 targetMapToOrigRow(i) = LO_INVALID;
771 targetMapToImportRow(i) = I_LID;
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);
784 template<
class Scalar,
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)
795 using Teuchos::Array;
796 using Teuchos::ArrayRCP;
797 using Teuchos::ArrayView;
802 typedef LocalOrdinal LO;
803 typedef GlobalOrdinal GO;
806 typedef Import<LO,GO,NO> import_type;
807 typedef Map<LO,GO,NO> map_type;
810 typedef typename map_type::local_map_type local_map_type;
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;
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")))));
823 LO LO_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
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();
838 lo_view_t Bcol2Ccol(Kokkos::ViewAllocateWithoutInitializing(
"Bcol2Ccol"),Pview.colMap->getNodeNumElements()), Icol2Ccol;
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));
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");
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));
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;
867 LO I_LID = Irowmap_local.getLocalElement(aidx);
868 targetMapToOrigRow(i) = LO_INVALID;
869 targetMapToImportRow(i) = I_LID;
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);
886 template<
class Scalar,
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,
895 const LocalOrdinalViewType & Acol2PIrow,
896 const LocalOrdinalViewType & Pcol2Accol,
897 const LocalOrdinalViewType & PIcol2Accol,
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"))));
908 using Teuchos::Array;
909 using Teuchos::ArrayRCP;
910 using Teuchos::ArrayView;
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;
923 typedef LocalOrdinal LO;
924 typedef GlobalOrdinal GO;
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();
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();
938 const KCRS & Amat = Aview.origMatrix->getLocalMatrix();
939 const KCRS & Pmat = Pview.origMatrix->getLocalMatrix();
940 const KCRS & Rmat = Rview.origMatrix->getLocalMatrix();
942 c_lno_view_t Arowptr = Amat.graph.row_map, Prowptr = Pmat.graph.row_map, Rrowptr = Rmat.graph.row_map;
943 const lno_nnz_view_t Acolind = Amat.graph.entries, Pcolind = Pmat.graph.entries , Rcolind = Rmat.graph.entries;
944 const scalar_view_t Avals = Amat.values, Pvals = Pmat.values, Rvals = Rmat.values;
946 c_lno_view_t Irowptr;
947 lno_nnz_view_t Icolind;
949 if(!Pview.importMatrix.is_null()) {
950 Irowptr = Pview.importMatrix->getLocalMatrix().graph.row_map;
951 Icolind = Pview.importMatrix->getLocalMatrix().graph.entries;
952 Ivals = Pview.importMatrix->getLocalMatrix().values;
953 p_max_nnz_per_row = std::max(p_max_nnz_per_row,Pview.importMatrix->getNodeMaxNumRowEntries());
956 #ifdef HAVE_TPETRA_MMM_TIMINGS
957 RCP<TimeMonitor> MM2 = rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string(
"RAP Newmatrix SerialCore - Compare"))));
967 size_t CSR_alloc = std::max(C_estimate_nnz(*Aview.origMatrix, *Pview.origMatrix), n);
968 lno_view_t Crowptr(Kokkos::ViewAllocateWithoutInitializing(
"Crowptr"),m+1);
969 lno_nnz_view_t Ccolind(Kokkos::ViewAllocateWithoutInitializing(
"Ccolind"),CSR_alloc);
970 scalar_view_t Cvals(Kokkos::ViewAllocateWithoutInitializing(
"Cvals"),CSR_alloc);
980 const size_t INVALID = Teuchos::OrdinalTraits<size_t>::invalid();
981 Array<size_t> ac_status(n, ST_INVALID);
991 size_t nnz = 0, nnz_old = 0;
992 for (
size_t i = 0; i < m; i++) {
998 for (
size_t kk = Rrowptr[i]; kk < Rrowptr[i+1]; kk++) {
1000 const SC Rik = Rvals[kk];
1004 for (
size_t ll = Arowptr[k]; ll < Arowptr[k+1]; ll++) {
1006 const SC Akl = Avals[ll];
1011 if (Acol2Prow[l] != LO_INVALID) {
1018 size_t Pl = Teuchos::as<size_t>(Acol2Prow[l]);
1021 for (
size_t jj = Prowptr[Pl]; jj < Prowptr[Pl+1]; jj++) {
1023 LO Acj = Pcol2Accol[j];
1026 if (ac_status[Acj] == INVALID || ac_status[Acj] < nnz_old) {
1027 #ifdef HAVE_TPETRA_DEBUG
1029 TEUCHOS_TEST_FOR_EXCEPTION(nnz >= Teuchos::as<size_t>(Ccolind.size()),
1031 label <<
" ERROR, not enough memory allocated for matrix product. Allocated: " << Ccolind.extent(0) << std::endl);
1034 ac_status[Acj] = nnz;
1036 Cvals[nnz] = Rik*Akl*Plj;
1039 Cvals[ac_status[Acj]] += Rik*Akl*Plj;
1049 size_t Il = Teuchos::as<size_t>(Acol2PIrow[l]);
1050 for (
size_t jj = Irowptr[Il]; jj < Irowptr[Il+1]; jj++) {
1052 LO Acj = PIcol2Accol[j];
1055 if (ac_status[Acj] == INVALID || ac_status[Acj] < nnz_old) {
1056 #ifdef HAVE_TPETRA_DEBUG
1058 TEUCHOS_TEST_FOR_EXCEPTION(nnz >= Teuchos::as<size_t>(Ccolind.size()),
1060 label <<
" ERROR, not enough memory allocated for matrix product. Allocated: " << Ccolind.extent(0) << std::endl);
1063 ac_status[Acj] = nnz;
1065 Cvals[nnz] = Rik*Akl*Plj;
1068 Cvals[ac_status[Acj]] += Rik*Akl*Plj;
1075 if (nnz + n > CSR_alloc) {
1077 Kokkos::resize(Ccolind,CSR_alloc);
1078 Kokkos::resize(Cvals,CSR_alloc);
1086 Kokkos::resize(Ccolind,nnz);
1087 Kokkos::resize(Cvals,nnz);
1089 #ifdef HAVE_TPETRA_MMM_TIMINGS
1090 MM = rcp(
new TimeMonitor (*TimeMonitor::getNewTimer(prefix_mmm + std::string(
"RAP Newmatrix Final Sort"))));
1094 if (params.is_null() || params->get(
"sort entries",
true))
1095 Import_Util::sortCrsEntries(Crowptr,Ccolind, Cvals);
1096 Ac.setAllValues(Crowptr, Ccolind, Cvals);
1098 #ifdef HAVE_TPETRA_MMM_TIMINGS
1099 MM = rcp(
new TimeMonitor (*TimeMonitor::getNewTimer(prefix_mmm + std::string(
"RAP Newmatrix ESFC"))));
1110 RCP<Teuchos::ParameterList> labelList = rcp(
new Teuchos::ParameterList);
1111 labelList->set(
"Timer Label",label);
1112 if(!params.is_null()) labelList->set(
"compute global constants",params->get(
"compute global constants",
true));
1113 RCP<const Export<LO,GO,NO> > dummyExport;
1114 Ac.expertStaticFillComplete(Pview.origMatrix->getDomainMap(),
1115 Rview.origMatrix->getRangeMap(),
1125 template<
class Scalar,
1127 class GlobalOrdinal,
1129 class LocalOrdinalViewType>
1130 void KernelWrappers3<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalOrdinalViewType>::mult_R_A_P_reuse_kernel_wrapper(CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Rview,
1131 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Aview,
1132 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Pview,
1133 const LocalOrdinalViewType & Acol2Prow,
1134 const LocalOrdinalViewType & Acol2PIrow,
1135 const LocalOrdinalViewType & Pcol2Accol,
1136 const LocalOrdinalViewType & PIcol2Accol,
1137 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Ac,
1138 Teuchos::RCP<
const Import<LocalOrdinal,GlobalOrdinal,Node> > Acimport,
1139 const std::string& label,
1140 const Teuchos::RCP<Teuchos::ParameterList>& params) {
1141 #ifdef HAVE_TPETRA_MMM_TIMINGS
1142 std::string prefix_mmm = std::string(
"TpetraExt ") + label + std::string(
": ");
1143 using Teuchos::TimeMonitor;
1144 Teuchos::RCP<Teuchos::TimeMonitor> MM = rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string(
"RAP Reuse SerialCore"))));
1147 using Teuchos::Array;
1148 using Teuchos::ArrayRCP;
1149 using Teuchos::ArrayView;
1155 typedef typename KCRS::StaticCrsGraphType graph_t;
1156 typedef typename graph_t::row_map_type::const_type c_lno_view_t;
1157 typedef typename graph_t::entries_type::non_const_type lno_nnz_view_t;
1158 typedef typename KCRS::values_type::non_const_type scalar_view_t;
1161 typedef LocalOrdinal LO;
1162 typedef GlobalOrdinal GO;
1164 typedef Map<LO,GO,NO> map_type;
1165 const size_t ST_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
1166 const LO LO_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
1167 const SC SC_ZERO = Teuchos::ScalarTraits<Scalar>::zero();
1170 RCP<const map_type> Accolmap = Ac.getColMap();
1171 size_t m = Rview.origMatrix->getNodeNumRows();
1172 size_t n = Accolmap->getNodeNumElements();
1173 size_t p_max_nnz_per_row = Pview.origMatrix->getNodeMaxNumRowEntries();
1176 const KCRS & Amat = Aview.origMatrix->getLocalMatrix();
1177 const KCRS & Pmat = Pview.origMatrix->getLocalMatrix();
1178 const KCRS & Rmat = Rview.origMatrix->getLocalMatrix();
1179 const KCRS & Cmat = Ac.getLocalMatrix();
1181 c_lno_view_t Arowptr = Amat.graph.row_map, Prowptr = Pmat.graph.row_map, Rrowptr = Rmat.graph.row_map, Crowptr = Cmat.graph.row_map;
1182 const lno_nnz_view_t Acolind = Amat.graph.entries, Pcolind = Pmat.graph.entries , Rcolind = Rmat.graph.entries, Ccolind = Cmat.graph.entries;
1183 const scalar_view_t Avals = Amat.values, Pvals = Pmat.values, Rvals = Rmat.values;
1184 scalar_view_t Cvals = Cmat.values;
1186 c_lno_view_t Irowptr;
1187 lno_nnz_view_t Icolind;
1188 scalar_view_t Ivals;
1189 if(!Pview.importMatrix.is_null()) {
1190 Irowptr = Pview.importMatrix->getLocalMatrix().graph.row_map;
1191 Icolind = Pview.importMatrix->getLocalMatrix().graph.entries;
1192 Ivals = Pview.importMatrix->getLocalMatrix().values;
1193 p_max_nnz_per_row = std::max(p_max_nnz_per_row,Pview.importMatrix->getNodeMaxNumRowEntries());
1196 #ifdef HAVE_TPETRA_MMM_TIMINGS
1197 RCP<TimeMonitor> MM2 = rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string(
"RAP Reuse SerialCore - Compare"))));
1208 Array<size_t> ac_status(n, ST_INVALID);
1218 size_t OLD_ip = 0, CSR_ip = 0;
1219 for (
size_t i = 0; i < m; i++) {
1222 OLD_ip = Crowptr[i];
1223 CSR_ip = Crowptr[i+1];
1224 for (
size_t k = OLD_ip; k < CSR_ip; k++) {
1225 ac_status[Ccolind[k]] = k;
1232 for (
size_t kk = Rrowptr[i]; kk < Rrowptr[i+1]; kk++) {
1234 const SC Rik = Rvals[kk];
1238 for (
size_t ll = Arowptr[k]; ll < Arowptr[k+1]; ll++) {
1240 const SC Akl = Avals[ll];
1245 if (Acol2Prow[l] != LO_INVALID) {
1252 size_t Pl = Teuchos::as<size_t>(Acol2Prow[l]);
1255 for (
size_t jj = Prowptr[Pl]; jj < Prowptr[Pl+1]; jj++) {
1257 LO Cij = Pcol2Accol[j];
1260 TEUCHOS_TEST_FOR_EXCEPTION(ac_status[Cij] < OLD_ip || ac_status[Cij] >= CSR_ip,
1261 std::runtime_error,
"Trying to insert a new entry (" << i <<
"," << Cij <<
") into a static graph " <<
1262 "(c_status = " << ac_status[Cij] <<
" of [" << OLD_ip <<
"," << CSR_ip <<
"))");
1264 Cvals[ac_status[Cij]] += Rik*Akl*Plj;
1273 size_t Il = Teuchos::as<size_t>(Acol2PIrow[l]);
1274 for (
size_t jj = Irowptr[Il]; jj < Irowptr[Il+1]; jj++) {
1276 LO Cij = PIcol2Accol[j];
1279 TEUCHOS_TEST_FOR_EXCEPTION(ac_status[Cij] < OLD_ip || ac_status[Cij] >= CSR_ip,
1280 std::runtime_error,
"Trying to insert a new entry (" << i <<
"," << Cij <<
") into a static graph " <<
1281 "(c_status = " << ac_status[Cij] <<
" of [" << OLD_ip <<
"," << CSR_ip <<
"))");
1283 Cvals[ac_status[Cij]] += Rik*Akl*Plj;
1290 #ifdef HAVE_TPETRA_MMM_TIMINGS
1291 auto MM3 = rcp(
new TimeMonitor (*TimeMonitor::getNewTimer(prefix_mmm + std::string(
"RAP Reuse ESFC"))));
1294 Ac.fillComplete(Ac.getDomainMap(), Ac.getRangeMap());
1302 template<
class Scalar,
1304 class GlobalOrdinal,
1306 class LocalOrdinalViewType>
1307 void KernelWrappers3<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalOrdinalViewType>::mult_PT_A_P_newmatrix_kernel_wrapper(CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Aview,
1308 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Pview,
1309 const LocalOrdinalViewType & Acol2Prow,
1310 const LocalOrdinalViewType & Acol2PIrow,
1311 const LocalOrdinalViewType & Pcol2Accol,
1312 const LocalOrdinalViewType & PIcol2Accol,
1313 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Ac,
1314 Teuchos::RCP<
const Import<LocalOrdinal,GlobalOrdinal,Node> > Acimport,
1315 const std::string& label,
1316 const Teuchos::RCP<Teuchos::ParameterList>& params) {
1317 #ifdef HAVE_TPETRA_MMM_TIMINGS
1318 std::string prefix_mmm = std::string(
"TpetraExt ") + label + std::string(
": ");
1319 using Teuchos::TimeMonitor;
1320 Teuchos::TimeMonitor MM(*TimeMonitor::getNewTimer(prefix_mmm + std::string(
"PTAP local transpose")));
1324 typedef RowMatrixTransposer<Scalar,LocalOrdinal,GlobalOrdinal, Node> transposer_type;
1325 transposer_type transposer (Pview.origMatrix,label+std::string(
"XP: "));
1327 using Teuchos::ParameterList;
1329 RCP<ParameterList> transposeParams (
new ParameterList);
1330 transposeParams->set (
"sort",
false);
1332 if (! params.is_null ()) {
1333 transposeParams->set (
"compute global constants",
1334 params->get (
"compute global constants: temporaries",
1337 RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Ptrans =
1338 transposer.createTransposeLocal (transposeParams);
1339 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node> Rview;
1340 Rview.origMatrix = Ptrans;
1342 mult_R_A_P_newmatrix_kernel_wrapper(Rview,Aview,Pview,Acol2Prow,Acol2PIrow,Pcol2Accol,PIcol2Accol,Ac,Acimport,label,params);
1348 template<
class Scalar,
1350 class GlobalOrdinal,
1352 class LocalOrdinalViewType>
1353 void KernelWrappers3<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalOrdinalViewType>::mult_PT_A_P_reuse_kernel_wrapper(CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Aview,
1354 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Pview,
1355 const LocalOrdinalViewType & Acol2Prow,
1356 const LocalOrdinalViewType & Acol2PIrow,
1357 const LocalOrdinalViewType & Pcol2Accol,
1358 const LocalOrdinalViewType & PIcol2Accol,
1359 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Ac,
1360 Teuchos::RCP<
const Import<LocalOrdinal,GlobalOrdinal,Node> > Acimport,
1361 const std::string& label,
1362 const Teuchos::RCP<Teuchos::ParameterList>& params) {
1363 #ifdef HAVE_TPETRA_MMM_TIMINGS
1364 std::string prefix_mmm = std::string(
"TpetraExt ") + label + std::string(
": ");
1365 using Teuchos::TimeMonitor;
1366 Teuchos::TimeMonitor MM(*TimeMonitor::getNewTimer(prefix_mmm + std::string(
"PTAP local transpose")));
1370 typedef RowMatrixTransposer<Scalar,LocalOrdinal,GlobalOrdinal, Node> transposer_type;
1371 transposer_type transposer (Pview.origMatrix,label+std::string(
"XP: "));
1373 using Teuchos::ParameterList;
1375 RCP<ParameterList> transposeParams (
new ParameterList);
1376 transposeParams->set (
"sort",
false);
1378 if (! params.is_null ()) {
1379 transposeParams->set (
"compute global constants",
1380 params->get (
"compute global constants: temporaries",
1383 RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Ptrans =
1384 transposer.createTransposeLocal (transposeParams);
1385 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node> Rview;
1386 Rview.origMatrix = Ptrans;
1388 mult_R_A_P_reuse_kernel_wrapper(Rview,Aview,Pview,Acol2Prow,Acol2PIrow,Pcol2Accol,PIcol2Accol,Ac,Acimport,label,params);
1396 template<
class Scalar,
1398 class GlobalOrdinal,
1400 void KernelWrappers3MMM<Scalar,LocalOrdinal,GlobalOrdinal,Node>::mult_PT_A_P_newmatrix_kernel_wrapper_2pass(CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Aview,
1401 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Pview,
1402 const Teuchos::Array<LocalOrdinal> & Acol2PRow,
1403 const Teuchos::Array<LocalOrdinal> & Acol2PRowImport,
1404 const Teuchos::Array<LocalOrdinal> & Pcol2Accol,
1405 const Teuchos::Array<LocalOrdinal> & PIcol2Accol,
1406 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Ac,
1407 Teuchos::RCP<
const Import<LocalOrdinal,GlobalOrdinal,Node> > Acimport,
1408 const std::string& label,
1409 const Teuchos::RCP<Teuchos::ParameterList>& params) {
1410 #ifdef HAVE_TPETRA_MMM_TIMINGS
1411 std::string prefix_mmm = std::string(
"TpetraExt ") + label + std::string(
": ");
1412 using Teuchos::TimeMonitor;
1413 Teuchos::RCP<Teuchos::TimeMonitor> MM = rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix_mmm + std::string(
"PTAP Newmatrix SerialCore"))));
1416 using Teuchos::Array;
1417 using Teuchos::ArrayRCP;
1418 using Teuchos::ArrayView;
1423 typedef LocalOrdinal LO;
1424 typedef GlobalOrdinal GO;
1426 typedef RowMatrixTransposer<SC,LO,GO,NO> transposer_type;
1427 const LO LO_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
1428 const SC SC_ZERO = Teuchos::ScalarTraits<Scalar>::zero();
1433 size_t n = Ac.getRowMap()->getNodeNumElements();
1434 LO maxAccol = Ac.getColMap()->getMaxLocalIndex();
1437 ArrayRCP<const size_t> Arowptr_RCP, Prowptr_RCP, Irowptr_RCP;
1438 ArrayRCP<size_t> Acrowptr_RCP;
1439 ArrayRCP<const LO> Acolind_RCP, Pcolind_RCP, Icolind_RCP;
1440 ArrayRCP<LO> Accolind_RCP;
1441 ArrayRCP<const Scalar> Avals_RCP, Pvals_RCP, Ivals_RCP;
1442 ArrayRCP<SC> Acvals_RCP;
1449 Aview.origMatrix->getAllValues(Arowptr_RCP, Acolind_RCP, Avals_RCP);
1450 Pview.origMatrix->getAllValues(Prowptr_RCP, Pcolind_RCP, Pvals_RCP);
1452 if (!Pview.importMatrix.is_null())
1453 Pview.importMatrix->getAllValues(Irowptr_RCP, Icolind_RCP, Ivals_RCP);
1460 ArrayView<const size_t> Arowptr, Prowptr, Irowptr;
1461 ArrayView<const LO> Acolind, Pcolind, Icolind;
1462 ArrayView<const SC> Avals, Pvals, Ivals;
1463 ArrayView<size_t> Acrowptr;
1464 ArrayView<LO> Accolind;
1465 ArrayView<SC> Acvals;
1466 Arowptr = Arowptr_RCP(); Acolind = Acolind_RCP(); Avals = Avals_RCP();
1467 Prowptr = Prowptr_RCP(); Pcolind = Pcolind_RCP(); Pvals = Pvals_RCP();
1468 if (!Pview.importMatrix.is_null()) {
1469 Irowptr = Irowptr_RCP(); Icolind = Icolind_RCP(); Ivals = Ivals_RCP();
1481 #ifdef HAVE_TPETRA_MMM_TIMINGS
1482 MM = rcp(
new TimeMonitor (*TimeMonitor::getNewTimer(prefix_mmm + std::string(
"PTAP local transpose"))));
1489 ArrayRCP<const size_t> Rrowptr_RCP;
1490 ArrayRCP<const LO> Rcolind_RCP;
1491 ArrayRCP<const Scalar> Rvals_RCP;
1492 ArrayView<const size_t> Rrowptr;
1493 ArrayView<const LO> Rcolind;
1494 ArrayView<const SC> Rvals;
1496 transposer_type transposer (Pview.origMatrix,label+std::string(
"XP: "));
1498 using Teuchos::ParameterList;
1499 RCP<ParameterList> transposeParams (
new ParameterList);
1500 transposeParams->set (
"sort",
false);
1501 if (! params.is_null ()) {
1502 transposeParams->set (
"compute global constants",
1503 params->get (
"compute global constants: temporaries",
1506 RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Ptrans =
1507 transposer.createTransposeLocal (transposeParams);
1509 Ptrans->getAllValues(Rrowptr_RCP, Rcolind_RCP, Rvals_RCP);
1510 Rrowptr = Rrowptr_RCP();
1511 Rcolind = Rcolind_RCP();
1512 Rvals = Rvals_RCP();
1517 #ifdef HAVE_TPETRA_MMM_TIMINGS
1518 MM = rcp(
new TimeMonitor (*TimeMonitor::getNewTimer(prefix_mmm + std::string(
"PTAP graph"))));
1521 const size_t ST_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
1522 Array<size_t> ac_status(maxAccol + 1, ST_INVALID);
1524 size_t nnz_alloc = std::max(Ac_estimate_nnz(*Aview.origMatrix, *Pview.origMatrix), n);
1525 size_t nnzPerRowA = 100;
1526 if (Aview.origMatrix->getNodeNumEntries() > 0)
1527 nnzPerRowA = Aview.origMatrix->getNodeNumEntries()/Aview.origMatrix->getNodeNumRows();
1528 Acrowptr_RCP.resize(n+1);
1529 Acrowptr = Acrowptr_RCP();
1530 Accolind_RCP.resize(nnz_alloc);
1531 Accolind = Accolind_RCP();
1533 size_t nnz = 0, nnz_old = 0;
1534 for (
size_t i = 0; i < n; i++) {
1540 for (
size_t kk = Rrowptr[i]; kk < Rrowptr[i+1]; kk++) {
1543 for (
size_t ll = Arowptr[k]; ll < Arowptr[k+1]; ll++) {
1546 if (Acol2PRow[l] != LO_INVALID) {
1553 size_t Pl = Teuchos::as<size_t>(Acol2PRow[l]);
1556 for (
size_t jj = Prowptr[Pl]; jj < Prowptr[Pl+1]; jj++) {
1558 LO Acj = Pcol2Accol[j];
1560 if (ac_status[Acj] == ST_INVALID || ac_status[Acj] < nnz_old) {
1562 ac_status[Acj] = nnz;
1563 Accolind[nnz] = Acj;
1574 size_t Il = Teuchos::as<size_t>(Acol2PRowImport[l]);
1575 for (
size_t jj = Irowptr[Il]; jj < Irowptr[Il+1]; jj++) {
1577 LO Acj = PIcol2Accol[j];
1579 if (ac_status[Acj] == ST_INVALID || ac_status[Acj] < nnz_old){
1581 ac_status[Acj] = nnz;
1582 Accolind[nnz] = Acj;
1592 if (nnz + std::max(5*nnzPerRowA, n) > nnz_alloc) {
1594 nnz_alloc = std::max(nnz_alloc, nnz + std::max(5*nnzPerRowA, n));
1595 Accolind_RCP.resize(nnz_alloc); Accolind = Accolind_RCP();
1596 Acvals_RCP.resize(nnz_alloc); Acvals = Acvals_RCP();
1603 Accolind_RCP.resize(nnz);
1604 Accolind = Accolind_RCP();
1607 Acvals_RCP.resize(nnz, SC_ZERO);
1608 Acvals = Acvals_RCP();
1615 #ifdef HAVE_TPETRA_MMM_TIMINGS
1616 MM = rcp(
new TimeMonitor (*TimeMonitor::getNewTimer(prefix_mmm + std::string(
"PTAP Newmatrix Fill Matrix"))));
1620 for (
size_t k = 0; k < n; k++) {
1621 for (
size_t ii = Prowptr[k]; ii < Prowptr[k+1]; ii++) {
1623 const SC Pki = Pvals[ii];
1624 for (
size_t ll = Arowptr[k]; ll < Arowptr[k+1]; ll++) {
1626 const SC Akl = Avals[ll];
1629 if (Acol2PRow[l] != LO_INVALID) {
1636 size_t Pl = Teuchos::as<size_t>(Acol2PRow[l]);
1637 for (
size_t jj = Prowptr[Pl]; jj < Prowptr[Pl+1]; jj++) {
1639 LO Acj = Pcol2Accol[j];
1641 for (pp = Acrowptr[i]; pp < Acrowptr[i+1]; pp++)
1642 if (Accolind[pp] == Acj)
1646 Acvals[pp] += Pki*Akl*Pvals[jj];
1655 size_t Il = Teuchos::as<size_t>(Acol2PRowImport[l]);
1656 for (
size_t jj = Irowptr[Il]; jj < Irowptr[Il+1]; jj++) {
1658 LO Acj = PIcol2Accol[j];
1660 for (pp = Acrowptr[i]; pp < Acrowptr[i+1]; pp++)
1661 if (Accolind[pp] == Acj)
1665 Acvals[pp] += Pki*Akl*Ivals[jj];
1673 #ifdef HAVE_TPETRA_MMM_TIMINGS
1674 MM = rcp(
new TimeMonitor (*TimeMonitor::getNewTimer(prefix_mmm + std::string(
"PTAP sort"))));
1681 Import_Util::sortCrsEntries(Acrowptr_RCP(), Accolind_RCP(), Acvals_RCP());
1684 Ac.setAllValues(Acrowptr_RCP, Accolind_RCP, Acvals_RCP);
1686 #ifdef HAVE_TPETRA_MMM_TIMINGS
1687 MM = rcp(
new TimeMonitor (*TimeMonitor::getNewTimer(prefix_mmm + std::string(
"PTAP Newmatrix ESFC"))));
1698 RCP<Teuchos::ParameterList> labelList = rcp(
new Teuchos::ParameterList);
1699 labelList->set(
"Timer Label",label);
1701 if(!params.is_null()) labelList->set(
"compute global constants",params->get(
"compute global constants",
true));
1702 RCP<const Export<LO,GO,NO> > dummyExport;
1703 Ac.expertStaticFillComplete(Pview.origMatrix->getDomainMap(),
1704 Pview.origMatrix->getDomainMap(),
1706 dummyExport, labelList);
1720 #define TPETRA_TRIPLEMATRIXMULTIPLY_INSTANT(SCALAR,LO,GO,NODE) \
1723 void TripleMatrixMultiply::MultiplyRAP( \
1724 const CrsMatrix< SCALAR , LO , GO , NODE >& R, \
1726 const CrsMatrix< SCALAR , LO , GO , NODE >& A, \
1728 const CrsMatrix< SCALAR , LO , GO , NODE >& P, \
1730 CrsMatrix< SCALAR , LO , GO , NODE >& Ac, \
1731 bool call_FillComplete_on_result, \
1732 const std::string & label, \
1733 const Teuchos::RCP<Teuchos::ParameterList>& params); \
1737 #endif // TPETRA_TRIPLEMATRIXMULTIPLY_DEF_HPP
Communication plan for data redistribution from a uniquely-owned to a (possibly) multiply-owned distr...
static int TAFC_OptimizationCoreCount()
The core count above which Tpetra::CrsMatrix::transferAndFillComplete will attempt to do advanced nei...
Sparse matrix that presents a row-oriented interface that lets users read or modify entries...
global_size_t getGlobalNumRows() const override
Number of global elements in the row map of this matrix.
bool isFillActive() const
Whether the matrix is not fill complete.
size_t global_size_t
Global size_t object.
bool isFillComplete() const override
Whether the matrix is fill complete.
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const override
The communicator over which the matrix is distributed.
Construct and (optionally) redistribute the explicitly stored transpose of a CrsMatrix.
Communication plan for data redistribution from a (possibly) multiply-owned to a uniquely-owned distr...
Teuchos::RCP< const RowGraph< LocalOrdinal, GlobalOrdinal, Node > > getGraph() const override
This matrix's graph, as a RowGraph.
Utility functions for packing and unpacking sparse matrix entries.
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 > ¶ms=Teuchos::null)
Sparse matrix-matrix multiply.
A parallel distribution of indices over processes.
Teuchos::RCP< const map_type > getDomainMap() const override
The domain Map of this matrix.
KokkosSparse::CrsMatrix< impl_scalar_type, local_ordinal_type, execution_space, void, typename local_graph_type::size_type > local_matrix_type
The specialization of Kokkos::CrsMatrix that represents the part of the sparse matrix on each MPI pro...
Internal functions and macros designed for use with Tpetra::Import and Tpetra::Export objects...
Stand-alone utility functions and macros.
Struct that holds views of the contents of a CrsMatrix.