Tpetra parallel linear algebra Version of the Day
Loading...
Searching...
No Matches
Tpetra_RowMatrixTransposer_def.hpp
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
42#ifndef TPETRA_ROWMATRIXTRANSPOSER_DEF_HPP
43#define TPETRA_ROWMATRIXTRANSPOSER_DEF_HPP
44
45#include "Tpetra_CrsMatrix.hpp"
46#include "Tpetra_BlockCrsMatrix.hpp"
47#include "Tpetra_Export.hpp"
50#include "Teuchos_ParameterList.hpp"
51#include "Teuchos_TimeMonitor.hpp"
52#include "KokkosSparse_Utils.hpp"
53#include "KokkosSparse_SortCrs.hpp"
54
55namespace Tpetra {
56
57template<class Scalar,
58 class LocalOrdinal,
59 class GlobalOrdinal,
60 class Node>
62RowMatrixTransposer (const Teuchos::RCP<const crs_matrix_type>& origMatrix,
63 const std::string& label)
64 : origMatrix_ (origMatrix), label_ (label)
65{}
66
67template<class Scalar,
68 class LocalOrdinal,
69 class GlobalOrdinal,
70 class Node>
71Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
73createTranspose (const Teuchos::RCP<Teuchos::ParameterList> &params)
74{
75 using Teuchos::RCP;
76 // Do the local transpose
77 RCP<crs_matrix_type> transMatrixWithSharedRows = createTransposeLocal (params);
78
79#ifdef HAVE_TPETRA_MMM_TIMINGS
80 const std::string prefix = std::string ("Tpetra ") + label_ + ": ";
81 using Teuchos::TimeMonitor;
82 TimeMonitor MM (*TimeMonitor::getNewTimer (prefix + "Transpose TAFC"));
83#endif
84
85 // If transMatrixWithSharedRows has an exporter, that's what we
86 // want. If it doesn't, the rows aren't actually shared, and we're
87 // done!
89 RCP<const export_type> exporter =
90 transMatrixWithSharedRows->getGraph ()->getExporter ();
91 if (exporter.is_null ()) {
92 return transMatrixWithSharedRows;
93 }
94 else {
95 Teuchos::ParameterList labelList;
96#ifdef HAVE_TPETRA_MMM_TIMINGS
97 labelList.set("Timer Label", label_);
98#endif
99 if(! params.is_null ()) {
100 const char paramName[] = "compute global constants";
101 labelList.set (paramName, params->get (paramName, true));
102 }
103 // Use the Export object to do a fused Export and fillComplete.
104 // This always sorts the local matrix after communication, so
105 // no need to set "sorted = false" in parameters.
106 return exportAndFillCompleteCrsMatrix<crs_matrix_type>
107 (transMatrixWithSharedRows, *exporter, Teuchos::null,
108 Teuchos::null, Teuchos::rcpFromRef (labelList));
109 }
110}
111
112
113// mfh 03 Feb 2013: In a definition outside the class like this, the
114// return value is considered outside the class scope (for things like
115// resolving typedefs), but the arguments are considered inside the
116// class scope.
117template<class Scalar,
118 class LocalOrdinal,
119 class GlobalOrdinal,
120 class Node>
121Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
123createTransposeLocal (const Teuchos::RCP<Teuchos::ParameterList>& params)
124{
125 using Teuchos::Array;
126 using Teuchos::ArrayRCP;
127 using Teuchos::ArrayView;
128 using Teuchos::RCP;
129 using Teuchos::rcp;
130 using Teuchos::rcp_dynamic_cast;
131 using LO = LocalOrdinal;
132 using GO = GlobalOrdinal;
133 using import_type = Tpetra::Import<LO, GO, Node>;
134 using export_type = Tpetra::Export<LO, GO, Node>;
135
136#ifdef HAVE_TPETRA_MMM_TIMINGS
137 std::string prefix = std::string("Tpetra ") + label_ + ": ";
138 using Teuchos::TimeMonitor;
139 TimeMonitor MM (*TimeMonitor::getNewTimer (prefix + "Transpose Local"));
140#endif
141
142 const bool sort = [&] () {
143 constexpr bool sortDefault = true; // see #4607 discussion
144 const char sortParamName[] = "sort";
145 return params.get () == nullptr ? sortDefault :
146 params->get (sortParamName, sortDefault);
147 } ();
148
149 const LO lclNumRows (origMatrix_->getLocalNumRows ());
150
151 RCP<const crs_matrix_type> crsMatrix =
152 rcp_dynamic_cast<const crs_matrix_type> (origMatrix_);
153 if (crsMatrix.is_null ()) {
154 auto rowMap = origMatrix_->getRowMap ();
155 if (rowMap->isOneToOne ()) {
156 Teuchos::Array<size_t> numEntPerRow (lclNumRows);
157 for (LO lclRow = 0; lclRow < lclNumRows; ++lclRow) {
158 numEntPerRow[lclRow] = origMatrix_->getNumEntriesInLocalRow (lclRow);
159 }
160 auto colMap = origMatrix_->getColMap ();
161
162 RCP<crs_matrix_type> crsMatrix_nc =
163 rcp (new crs_matrix_type (rowMap, colMap, numEntPerRow ()));
164
165 // When source & target Maps are same, Import just copies.
166 import_type imp (rowMap, rowMap);
167 crsMatrix_nc->doImport (*origMatrix_, imp, Tpetra::REPLACE);
168 crsMatrix_nc->fillComplete (origMatrix_->getDomainMap (),
169 origMatrix_->getRangeMap ());
170 crsMatrix = crsMatrix_nc;
171 }
172 else {
173 TEUCHOS_ASSERT( false ); // not implemented (it wasn't before)
174 }
175 }
176
177 using local_matrix_device_type = typename crs_matrix_type::local_matrix_device_type;
178
179 local_matrix_device_type lclMatrix = crsMatrix->getLocalMatrixDevice ();
180 local_matrix_device_type lclTransposeMatrix = KokkosSparse::Impl::transpose_matrix(lclMatrix);
181 if (sort)
182 KokkosSparse::sort_crs_matrix(lclTransposeMatrix);
183
184 // Prebuild the importers and exporters the no-communication way,
185 // flipping the importers and exporters around.
186 const auto origExport = origMatrix_->getGraph ()->getExporter ();
187 RCP<const import_type> myImport = origExport.is_null () ?
188 Teuchos::null : rcp (new import_type (*origExport));
189 const auto origImport = origMatrix_->getGraph ()->getImporter ();
190 RCP<const export_type> myExport = origImport.is_null () ?
191 Teuchos::null : rcp (new export_type (*origImport));
192
193 RCP<Teuchos::ParameterList> graphParams = Teuchos::null;
194 if(!sort) {
195 graphParams = rcp(new Teuchos::ParameterList);
196 graphParams->set("sorted", false);
197 }
198
199 return rcp (new crs_matrix_type (lclTransposeMatrix,
200 origMatrix_->getColMap (),
201 origMatrix_->getRowMap (),
202 origMatrix_->getRangeMap (),
203 origMatrix_->getDomainMap (),
204 myImport, myExport, graphParams));
205}
206
207/*************************************************************************/
208
209template<class Scalar,
210 class LocalOrdinal,
211 class GlobalOrdinal,
212 class Node>
214BlockCrsMatrixTransposer (const Teuchos::RCP<const bcrs_matrix_type>& origMatrix,
215 const std::string& label)
216 : origMatrix_ (origMatrix), label_ (label)
217{}
218
219template<class Scalar,
220 class LocalOrdinal,
221 class GlobalOrdinal,
222 class Node>
223Teuchos::RCP<BlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
225createTranspose (const Teuchos::RCP<Teuchos::ParameterList> &params)
226{
227 using Teuchos::RCP;
228 // Do the local transpose
229 RCP<bcrs_matrix_type> transMatrixWithSharedRows = createTransposeLocal (params);
230
231#ifdef HAVE_TPETRA_MMM_TIMINGS
232 const std::string prefix = std::string ("Tpetra ") + label_ + ": ";
233 using Teuchos::TimeMonitor;
234 TimeMonitor MM (*TimeMonitor::getNewTimer (prefix + "Transpose TAFC"));
235#endif
236
237 // If transMatrixWithSharedRows has an exporter, that's what we
238 // want. If it doesn't, the rows aren't actually shared, and we're
239 // done!
241 RCP<const export_type> exporter =
242 transMatrixWithSharedRows->getGraph ()->getExporter ();
243 if (exporter.is_null ()) {
244 return transMatrixWithSharedRows;
245 }
246 else {
247 Teuchos::ParameterList labelList;
248#ifdef HAVE_TPETRA_MMM_TIMINGS
249 labelList.set("Timer Label", label_);
250#endif
251 if(! params.is_null ()) {
252 const char paramName[] = "compute global constants";
253 labelList.set (paramName, params->get (paramName, true));
254 }
255 // Use the Export object to do a fused Export and fillComplete.
256 // This always sorts the local matrix after communication, so
257 // no need to set "sorted = false" in parameters.
258 return exportAndFillCompleteBlockCrsMatrix<bcrs_matrix_type>
259 (transMatrixWithSharedRows, *exporter);
260 }
261}
262
263
264// mfh 03 Feb 2013: In a definition outside the class like this, the
265// return value is considered outside the class scope (for things like
266// resolving typedefs), but the arguments are considered inside the
267// class scope.
268template<class Scalar,
269 class LocalOrdinal,
270 class GlobalOrdinal,
271 class Node>
272Teuchos::RCP<BlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
274createTransposeLocal (const Teuchos::RCP<Teuchos::ParameterList>& params)
275{
276 using Teuchos::Array;
277 using Teuchos::ArrayRCP;
278 using Teuchos::ArrayView;
279 using Teuchos::RCP;
280 using Teuchos::rcp;
281 using Teuchos::rcp_dynamic_cast;
282 using LO = LocalOrdinal;
283 using GO = GlobalOrdinal;
284 using import_type = Tpetra::Import<LO, GO, Node>;
285 using export_type = Tpetra::Export<LO, GO, Node>;
286 using crs_graph_type = typename bcrs_matrix_type::crs_graph_type;
287
288#ifdef HAVE_TPETRA_MMM_TIMINGS
289 std::string prefix = std::string("Tpetra ") + label_ + ": ";
290 using Teuchos::TimeMonitor;
291 TimeMonitor MM (*TimeMonitor::getNewTimer (prefix + "Transpose Local"));
292#endif
293
294 RCP<const bcrs_matrix_type> crsMatrix =
295 rcp_dynamic_cast<const bcrs_matrix_type> (origMatrix_);
296
297 if(crsMatrix.is_null())
298 TEUCHOS_ASSERT( false ); // not implemented
299
300 using local_matrix_device_type = typename bcrs_matrix_type::local_matrix_device_type;
301
302 typename local_matrix_device_type::values_type values ;
303 RCP<const crs_graph_type> graph;
304 {
305 local_matrix_device_type lclMatrix = crsMatrix->getLocalMatrixDevice ();
306
307 local_matrix_device_type lclTransposeMatrix = KokkosSparse::Impl::transpose_bsr_matrix(lclMatrix);
308
309 // BlockCrs requires that we sort stuff
310 KokkosSparse::sort_crs_matrix(lclTransposeMatrix);
311 values = lclTransposeMatrix.values;
312
313 // Prebuild the importers and exporters the no-communication way,
314 // flipping the importers and exporters around.
315 const auto origExport = origMatrix_->getGraph ()->getExporter ();
316 RCP<const import_type> myImport = origExport.is_null () ?
317 Teuchos::null : rcp (new import_type (*origExport));
318 const auto origImport = origMatrix_->getGraph ()->getImporter ();
319 RCP<const export_type> myExport = origImport.is_null () ?
320 Teuchos::null : rcp (new export_type (*origImport));
321
322 RCP<Teuchos::ParameterList> graphParams = Teuchos::null;
323
324 // Make the Transpose Graph
325 graph = rcp(new crs_graph_type(lclTransposeMatrix.graph,
326 origMatrix_->getColMap (),
327 origMatrix_->getRowMap (),
328 origMatrix_->getGraph()->getRangeMap (),
329 origMatrix_->getGraph()->getDomainMap (),
330 myImport,
331 myExport,
332 graphParams));
333 }
334 // Now make the matrix
335 return rcp (new bcrs_matrix_type (*graph,
336 values,
337 origMatrix_->getBlockSize()));
338}
339//
340
341
342//
343// Explicit instantiation macro
344//
345// Must be expanded from within the Tpetra namespace!
346//
347
348#define TPETRA_ROWMATRIXTRANSPOSER_INSTANT(SCALAR,LO,GO,NODE) \
349 template class RowMatrixTransposer< SCALAR, LO , GO , NODE >;\
350 template class BlockCrsMatrixTransposer< SCALAR, LO , GO , NODE >;
351
352} // namespace Tpetra
353
354#endif
Declare and define the functions Tpetra::Details::computeOffsetsFromCounts and Tpetra::computeOffsets...
Declaration and definition of functions for sorting "short" arrays of keys and corresponding values.
Teuchos::RCP< bcrs_matrix_type > createTransposeLocal(const Teuchos::RCP< Teuchos::ParameterList > &params=Teuchos::null)
Compute and return the transpose of the matrix given to the constructor.
BlockCrsMatrixTransposer(const Teuchos::RCP< const bcrs_matrix_type > &origMatrix, const std::string &label=std::string())
Constructor that takes the matrix to transpose.
Teuchos::RCP< bcrs_matrix_type > createTranspose(const Teuchos::RCP< Teuchos::ParameterList > &params=Teuchos::null)
Compute and return the transpose of the matrix given to the constructor.
Sparse matrix whose entries are small dense square blocks, all of the same dimensions.
Sparse matrix that presents a row-oriented interface that lets users read or modify entries.
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...
Communication plan for data redistribution from a (possibly) multiply-owned to a uniquely-owned distr...
Communication plan for data redistribution from a uniquely-owned to a (possibly) multiply-owned distr...
Teuchos::RCP< crs_matrix_type > createTranspose(const Teuchos::RCP< Teuchos::ParameterList > &params=Teuchos::null)
Compute and return the transpose of the matrix given to the constructor.
Teuchos::RCP< crs_matrix_type > createTransposeLocal(const Teuchos::RCP< Teuchos::ParameterList > &params=Teuchos::null)
Compute and return the transpose of the matrix given to the constructor.
RowMatrixTransposer(const Teuchos::RCP< const crs_matrix_type > &origMatrix, const std::string &label=std::string())
Constructor that takes the matrix to transpose.
Namespace Tpetra contains the class and methods constituting the Tpetra library.
void sort(View &view, const size_t &size)
Convenience wrapper for std::sort for host-accessible views.
@ REPLACE
Replace existing values with new values.