Tpetra parallel linear algebra  Version of the Day
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 
46 #include <Tpetra_CrsMatrix.hpp>
47 #include <Tpetra_Export.hpp>
48 #include <Tpetra_Import.hpp>
49 #include <Teuchos_TimeMonitor.hpp>
50 
51 namespace Tpetra {
52 
53 template<class Scalar,
54  class LocalOrdinal,
55  class GlobalOrdinal,
56  class Node>
58 RowMatrixTransposer (const Teuchos::RCP<const crs_matrix_type>& origMatrix,const std::string & label)
59  : origMatrix_(origMatrix), label_(label) {}
60 
61 template<class Scalar,
62  class LocalOrdinal,
63  class GlobalOrdinal,
64  class Node>
65 Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
68 {
69  using Teuchos::RCP;
70 #ifdef HAVE_TPETRA_MMM_TIMINGS
71  std::string prefix = std::string("Tpetra ")+ label_ + std::string(": ");
72  using Teuchos::TimeMonitor;
73  Teuchos::RCP<Teuchos::TimeMonitor> MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string("Transpose Local"))));
74 #endif
75  // Do the local transpose
76  RCP<crs_matrix_type> transMatrixWithSharedRows = createTransposeLocal ();
77 
78 #ifdef HAVE_TPETRA_MMM_TIMINGS
79  MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string("Transpose TAFC"))));
80 #endif
81 
82  // If transMatrixWithSharedRows has an exporter, that's what we
83  // want. If it doesn't, the rows aren't actually shared, and we're
84  // done!
85  RCP<const Export<LocalOrdinal,GlobalOrdinal,Node> > exporter =
86  transMatrixWithSharedRows->getGraph ()->getExporter ();
87  if (exporter.is_null ()) {
88  return transMatrixWithSharedRows;
89  }
90  else {
91  Teuchos::ParameterList labelList;
92 #ifdef HAVE_TPETRA_MMM_TIMINGS
93  labelList.set("Timer Label",label_);
94 #endif
95  // Use the Export object to do a fused Export and fillComplete.
96  return exportAndFillCompleteCrsMatrix<crs_matrix_type> (transMatrixWithSharedRows, *exporter,Teuchos::null,Teuchos::null,Teuchos::rcp(&labelList,false));
97  }
98 }
99 
100 
101 // mfh 03 Feb 2013: In a definition outside the class like this, the
102 // return value is considered outside the class scope (for things like
103 // resolving typedefs), but the arguments are considered inside the
104 // class scope.
105 template<class Scalar,
106  class LocalOrdinal,
107  class GlobalOrdinal,
108  class Node>
109 Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
112 {
113  using Teuchos::Array;
114  using Teuchos::ArrayRCP;
115  using Teuchos::ArrayView;
116  using Teuchos::RCP;
117  using Teuchos::rcp;
118  using Teuchos::rcp_dynamic_cast;
119  typedef LocalOrdinal LO;
120  typedef GlobalOrdinal GO;
121  typedef Tpetra::Import<LO, GO, Node> import_type;
122  typedef Tpetra::Export<LO, GO, Node> export_type;
123 
124  //
125  // This transpose is based upon the approach in EpetraExt.
126  //
127  size_t numLocalCols = origMatrix_->getNodeNumCols();
128  size_t numLocalRows = origMatrix_->getNodeNumRows();
129  size_t numLocalNnz = origMatrix_->getNodeNumEntries();
130 
131  // Determine how many nonzeros there are per row in the transpose.
132  Array<size_t> CurrentStart(numLocalCols,0);
133  ArrayView<const LO> localIndices;
134  ArrayView<const Scalar> localValues;
135  RCP<const crs_matrix_type> crsMatrix =
136  rcp_dynamic_cast<const crs_matrix_type> (origMatrix_);
137  if (crsMatrix == Teuchos::null) {
138  for (size_t i=0; i<numLocalRows; ++i) {
139  const size_t numEntriesInRow = origMatrix_->getNumEntriesInLocalRow(i);
140  origMatrix_->getLocalRowView(i, localIndices, localValues);
141  for (size_t j=0; j<numEntriesInRow; ++j) {
142  ++CurrentStart[ localIndices[j] ];
143  }
144  }
145  } else {
146  ArrayRCP<const size_t> origRowPtr_rcp;
147  ArrayRCP<const LO> origColInd_rcp;
148  ArrayRCP<const Scalar> origValues_rcp;
149  crsMatrix->getAllValues(origRowPtr_rcp, origColInd_rcp, origValues_rcp);
150  ArrayView<const LO> origColInd = origColInd_rcp();
151  for (LO j=0; j<origColInd.size(); ++j) {
152  ++CurrentStart[ origColInd[j] ];
153  }
154  }
155 
156  // create temporary row-major storage for the transposed matrix
157 
158  ArrayRCP<size_t> rowptr_rcp(numLocalCols+1);
159  ArrayRCP<LO> colind_rcp(numLocalNnz);
160  ArrayRCP<Scalar> values_rcp(numLocalNnz);
161 
162  // Since ArrayRCP's are slow...
163  ArrayView<size_t> TransRowptr = rowptr_rcp();
164  ArrayView<LO> TransColind = colind_rcp();
165  ArrayView<Scalar> TransValues = values_rcp();
166 
167  // Scansum the TransRowptr; reset CurrentStart
168  TransRowptr[0]=0;
169  for (size_t i=1; i<numLocalCols+1; ++i) TransRowptr[i] = CurrentStart[i-1] + TransRowptr[i-1];
170  for (size_t i=0; i<numLocalCols; ++i) CurrentStart[i] = TransRowptr[i];
171 
172  // populate the row-major storage so that the data for the transposed
173  // matrix is easy to access
174  if (crsMatrix == Teuchos::null) {
175  for (size_t i=0; i<numLocalRows; ++i) {
176  const size_t numEntriesInRow = origMatrix_->getNumEntriesInLocalRow (i);
177  origMatrix_->getLocalRowView(i, localIndices, localValues);
178 
179  for (size_t j=0; j<numEntriesInRow; ++j) {
180  size_t idx = CurrentStart[localIndices[j]];
181  TransColind[idx] = Teuchos::as<LO>(i);
182  TransValues[idx] = localValues[j];
183  ++CurrentStart[localIndices[j]];
184  }
185  } //for (size_t i=0; i<numLocalRows; ++i)
186  } else {
187  ArrayRCP<const size_t> origRowPtr_rcp;
188  ArrayRCP<const LO> origColInd_rcp;
189  ArrayRCP<const Scalar> origValues_rcp;
190  crsMatrix->getAllValues(origRowPtr_rcp, origColInd_rcp, origValues_rcp);
191  ArrayView<const size_t> origRowPtr = origRowPtr_rcp();
192  ArrayView<const LO> origColInd = origColInd_rcp();
193  ArrayView<const Scalar> origValues = origValues_rcp();
194  size_t k=0;
195  for (LO i=0; i<origRowPtr.size()-1; ++i) {
196  const LO rowIndex = Teuchos::as<LO>(i);
197  size_t rowStart = origRowPtr[i], rowEnd = origRowPtr[i+1];
198  for (size_t j = rowStart; j < rowEnd; ++j) {
199  size_t idx = CurrentStart[origColInd[k]];
200  TransColind[idx] = rowIndex;
201  TransValues[idx] = origValues[k];
202  ++CurrentStart[origColInd[k++]];
203  }
204  }
205  }
206 
207  //Allocate and populate temporary matrix with rows not uniquely owned
208  RCP<crs_matrix_type> transMatrixWithSharedRows =
209  rcp (new crs_matrix_type (origMatrix_->getColMap (),
210  origMatrix_->getRowMap (), 0));
211  transMatrixWithSharedRows->setAllValues (rowptr_rcp, colind_rcp, values_rcp);
212 
213  // Prebuild the importers and exporters the no-communication way,
214  // flipping the importers and exporters around.
215  RCP<const import_type> myImport;
216  RCP<const export_type> myExport;
217  if (! origMatrix_->getGraph ()->getImporter ().is_null ()) {
218  myExport = rcp (new export_type (*origMatrix_->getGraph ()->getImporter ()));
219  }
220  if (! origMatrix_->getGraph ()->getExporter ().is_null ()) {
221  myImport = rcp (new import_type (*origMatrix_->getGraph ()->getExporter ()));
222  }
223 
224  // Call ESFC & return
225  transMatrixWithSharedRows->expertStaticFillComplete (origMatrix_->getRangeMap (),
226  origMatrix_->getDomainMap (),
227  myImport, myExport);
228  return transMatrixWithSharedRows;
229 }
230 //
231 // Explicit instantiation macro
232 //
233 // Must be expanded from within the Tpetra namespace!
234 //
235 
236 #define TPETRA_ROWMATRIXTRANSPOSER_INSTANT(SCALAR,LO,GO,NODE) \
237  \
238  template class RowMatrixTransposer< SCALAR, LO , GO , NODE >;
239 
240 
241 }
242 
243 #endif
RowMatrixTransposer(const Teuchos::RCP< const crs_matrix_type > &origMatrix, const std::string &label=std::string())
Constructor that takes the matrix to transpose.
Communication plan for data redistribution from a uniquely-owned to a (possibly) multiply-owned distr...
Namespace Tpetra contains the class and methods constituting the Tpetra library.
Sparse matrix that presents a row-oriented interface that lets users read or modify entries...
Communication plan for data redistribution from a (possibly) multiply-owned to a uniquely-owned distr...
Teuchos::RCP< crs_matrix_type > createTranspose()
Compute and return the transpose of the matrix given to the constructor.
size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const
Returns the current number of entries on this node in the specified local row.
Teuchos::RCP< crs_matrix_type > createTransposeLocal()
Compute and return the transpose of the matrix given to the constructor.