Teko  Version of the Day
 All Classes Files Functions Variables Pages
Teko_BlockingTpetra.cpp
1 // @HEADER
2 // *****************************************************************************
3 // Teko: A package for block and physics based preconditioning
4 //
5 // Copyright 2010 NTESS and the Teko contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #include "Teko_BlockingTpetra.hpp"
11 #include "Teko_Utilities.hpp"
12 
13 #include "Tpetra_Vector.hpp"
14 #include "Tpetra_Map.hpp"
15 
16 using Teuchos::RCP;
17 using Teuchos::rcp;
18 
19 namespace Teko {
20 namespace TpetraHelpers {
21 namespace Blocking {
22 
36 const MapPair buildSubMap(const std::vector<GO>& gid, const Teuchos::Comm<int>& comm) {
37  using GST = Tpetra::global_size_t;
38  const GST invalid = Teuchos::OrdinalTraits<GST>::invalid();
39  Teuchos::RCP<Tpetra::Map<LO, GO, NT> > gidMap = rcp(
40  new Tpetra::Map<LO, GO, NT>(invalid, Teuchos::ArrayView<const GO>(gid), 0, rcpFromRef(comm)));
41  Teuchos::RCP<Tpetra::Map<LO, GO, NT> > contigMap =
42  rcp(new Tpetra::Map<LO, GO, NT>(invalid, gid.size(), 0, rcpFromRef(comm)));
43 
44  return std::make_pair(gidMap, contigMap);
45 }
46 
56 const ImExPair buildExportImport(const Tpetra::Map<LO, GO, NT>& baseMap, const MapPair& maps) {
57  return std::make_pair(rcp(new Tpetra::Import<LO, GO, NT>(rcpFromRef(baseMap), maps.first)),
58  rcp(new Tpetra::Export<LO, GO, NT>(maps.first, rcpFromRef(baseMap))));
59 }
60 
68 void buildSubVectors(const std::vector<MapPair>& maps,
69  std::vector<RCP<Tpetra::MultiVector<ST, LO, GO, NT> > >& vectors, int count) {
70  std::vector<MapPair>::const_iterator mapItr;
71 
72  // loop over all maps
73  for (mapItr = maps.begin(); mapItr != maps.end(); mapItr++) {
74  // add new elements to vectors
75  RCP<Tpetra::MultiVector<ST, LO, GO, NT> > mv =
76  rcp(new Tpetra::MultiVector<ST, LO, GO, NT>((*mapItr).second, count));
77  vectors.push_back(mv);
78  }
79 }
80 
88 void one2many(std::vector<RCP<Tpetra::MultiVector<ST, LO, GO, NT> > >& many,
89  const Tpetra::MultiVector<ST, LO, GO, NT>& single,
90  const std::vector<RCP<Tpetra::Import<LO, GO, NT> > >& subImport) {
91  // std::vector<RCP<Epetra_Vector> >::const_iterator vecItr;
92  std::vector<RCP<Tpetra::MultiVector<ST, LO, GO, NT> > >::const_iterator vecItr;
93  std::vector<RCP<Tpetra::Import<LO, GO, NT> > >::const_iterator impItr;
94 
95  // using Importers fill the sub vectors from the mama vector
96  for (vecItr = many.begin(), impItr = subImport.begin(); vecItr != many.end();
97  ++vecItr, ++impItr) {
98  // for ease of access to the destination
99  RCP<Tpetra::MultiVector<ST, LO, GO, NT> > destVec = *vecItr;
100 
101  // extract the map with global indicies from the current vector
102  const Tpetra::Map<LO, GO, NT>& globalMap = *(*impItr)->getTargetMap();
103 
104  // build the import vector as a view on the destination
105  GO lda = destVec->getStride();
106  GO destSize = destVec->getGlobalLength() * destVec->getNumVectors();
107  std::vector<ST> destArray(destSize);
108  Teuchos::ArrayView<ST> destVals(destArray);
109  destVec->get1dCopy(destVals, lda);
110  Tpetra::MultiVector<ST, LO, GO, NT> importVector(rcpFromRef(globalMap), destVals, lda,
111  destVec->getNumVectors());
112 
113  // perform the import
114  importVector.doImport(single, **impItr, Tpetra::INSERT);
115 
116  Tpetra::Import<LO, GO, NT> importer(destVec->getMap(), destVec->getMap());
117  importVector.replaceMap(destVec->getMap());
118  destVec->doExport(importVector, importer, Tpetra::INSERT);
119  }
120 }
121 
132 void many2one(Tpetra::MultiVector<ST, LO, GO, NT>& one,
133  const std::vector<RCP<const Tpetra::MultiVector<ST, LO, GO, NT> > >& many,
134  const std::vector<RCP<Tpetra::Export<LO, GO, NT> > >& subExport) {
135  // std::vector<RCP<const Epetra_Vector> >::const_iterator vecItr;
136  std::vector<RCP<const Tpetra::MultiVector<ST, LO, GO, NT> > >::const_iterator vecItr;
137  std::vector<RCP<Tpetra::Export<LO, GO, NT> > >::const_iterator expItr;
138 
139  // using Exporters fill the empty vector from the sub-vectors
140  for (vecItr = many.begin(), expItr = subExport.begin(); vecItr != many.end();
141  ++vecItr, ++expItr) {
142  // for ease of access to the source
143  RCP<const Tpetra::MultiVector<ST, LO, GO, NT> > srcVec = *vecItr;
144 
145  // extract the map with global indicies from the current vector
146  const Tpetra::Map<LO, GO, NT>& globalMap = *(*expItr)->getSourceMap();
147 
148  // build the export vector as a view of the destination
149  GO lda = srcVec->getStride();
150  GO srcSize = srcVec->getGlobalLength() * srcVec->getNumVectors();
151  std::vector<ST> srcArray(srcSize);
152  Teuchos::ArrayView<ST> srcVals(srcArray);
153  srcVec->get1dCopy(srcVals, lda);
154  Tpetra::MultiVector<ST, LO, GO, NT> exportVector(rcpFromRef(globalMap), srcVals, lda,
155  srcVec->getNumVectors());
156 
157  // perform the export
158  one.doExport(exportVector, **expItr, Tpetra::INSERT);
159  }
160 }
161 
167 RCP<Tpetra::Vector<GO, LO, GO, NT> > getSubBlockColumnGIDs(
168  const Tpetra::CrsMatrix<ST, LO, GO, NT>& A, const MapPair& mapPair) {
169  RCP<const Tpetra::Map<LO, GO, NT> > blkGIDMap = mapPair.first;
170  RCP<const Tpetra::Map<LO, GO, NT> > blkContigMap = mapPair.second;
171 
172  // fill index vector for rows
173  Tpetra::Vector<GO, LO, GO, NT> rIndex(A.getRowMap(), true);
174  for (size_t i = 0; i < A.getLocalNumRows(); i++) {
175  // LID is need to get to contiguous map
176  LO lid = blkGIDMap->getLocalElement(
177  A.getRowMap()->getGlobalElement(i)); // this LID makes me nervous
178  if (lid > -1)
179  rIndex.replaceLocalValue(i, blkContigMap->getGlobalElement(lid));
180  else
181  rIndex.replaceLocalValue(i, -1);
182  }
183 
184  // get relavant column indices
185 
186  Tpetra::Import<LO, GO, NT> import(A.getRowMap(), A.getColMap());
187  RCP<Tpetra::Vector<GO, LO, GO, NT> > cIndex =
188  rcp(new Tpetra::Vector<GO, LO, GO, NT>(A.getColMap(), true));
189  cIndex->doImport(rIndex, import, Tpetra::INSERT);
190 
191  return cIndex;
192 }
193 
194 // build a single subblock Epetra_CrsMatrix
195 RCP<Tpetra::CrsMatrix<ST, LO, GO, NT> > buildSubBlock(
196  int i, int j, const RCP<const Tpetra::CrsMatrix<ST, LO, GO, NT> >& A,
197  const std::vector<MapPair>& subMaps) {
198  // get the number of variables families
199  int numVarFamily = subMaps.size();
200 
201  TEUCHOS_ASSERT(i >= 0 && i < numVarFamily);
202  TEUCHOS_ASSERT(j >= 0 && j < numVarFamily);
203  const RCP<const Tpetra::Map<LO, GO, NT> > gRowMap = subMaps[i].first; // new GIDs
204  const RCP<const Tpetra::Map<LO, GO, NT> > rowMap = subMaps[i].second; // contiguous GIDs
205  const RCP<const Tpetra::Map<LO, GO, NT> > colMap = subMaps[j].second;
206 
207  const RCP<Tpetra::Vector<GO, LO, GO, NT> > plocal2ContigGIDs =
208  getSubBlockColumnGIDs(*A, subMaps[j]);
209  Tpetra::Vector<GO, LO, GO, NT>& local2ContigGIDs = *plocal2ContigGIDs;
210 
211  // get entry information
212  LO numMyRows = rowMap->getLocalNumElements();
213  LO maxNumEntries = A->getLocalMaxNumRowEntries();
214 
215  // for extraction
216  auto indices = typename Tpetra::CrsMatrix<ST, LO, GO, NT>::nonconst_local_inds_host_view_type(
217  Kokkos::ViewAllocateWithoutInitializing("rowIndices"), maxNumEntries);
218  auto values = typename Tpetra::CrsMatrix<ST, LO, GO, NT>::nonconst_values_host_view_type(
219  Kokkos::ViewAllocateWithoutInitializing("rowIndices"), maxNumEntries);
220 
221  // for insertion
222  std::vector<GO> colIndices(maxNumEntries);
223  std::vector<ST> colValues(maxNumEntries);
224 
225  // for counting
226  std::vector<size_t> nEntriesPerRow(numMyRows);
227 
228  const size_t invalid = Teuchos::OrdinalTraits<size_t>::invalid();
229 
230  // insert each row into subblock
231  // Count the number of entries per row in the new matrix
232  for (LO localRow = 0; localRow < numMyRows; localRow++) {
233  size_t numEntries = invalid;
234  GO globalRow = gRowMap->getGlobalElement(localRow);
235  LO lid = A->getRowMap()->getLocalElement(globalRow);
236  TEUCHOS_ASSERT(lid > -1);
237 
238  A->getLocalRowCopy(lid, indices, values, numEntries);
239 
240  LO numOwnedCols = 0;
241  for (size_t localCol = 0; localCol < numEntries; localCol++) {
242  TEUCHOS_ASSERT(indices(localCol) > -1);
243 
244  // if global id is not owned by this column
245 
246  GO gid = local2ContigGIDs.getData()[indices[localCol]];
247  if (gid == -1) continue; // in contiguous row
248  numOwnedCols++;
249  }
250  nEntriesPerRow[localRow] = numOwnedCols;
251  }
252 
253  RCP<Tpetra::CrsMatrix<ST, LO, GO, NT> > mat = rcp(
254  new Tpetra::CrsMatrix<ST, LO, GO, NT>(rowMap, Teuchos::ArrayView<size_t>(nEntriesPerRow)));
255 
256  // insert each row into subblock
257  // let FillComplete handle column distribution
258  for (LO localRow = 0; localRow < numMyRows; localRow++) {
259  size_t numEntries = invalid;
260  GO globalRow = gRowMap->getGlobalElement(localRow);
261  LO lid = A->getRowMap()->getLocalElement(globalRow);
262  GO contigRow = rowMap->getGlobalElement(localRow);
263  TEUCHOS_ASSERT(lid > -1);
264 
265  A->getLocalRowCopy(lid, indices, values, numEntries);
266 
267  LO numOwnedCols = 0;
268  for (size_t localCol = 0; localCol < numEntries; localCol++) {
269  TEUCHOS_ASSERT(indices(localCol) > -1);
270 
271  // if global id is not owned by this column
272 
273  GO gid = local2ContigGIDs.getData()[indices(localCol)];
274  if (gid == -1) continue; // in contiguous row
275 
276  colIndices[numOwnedCols] = gid;
277  colValues[numOwnedCols] = values(localCol);
278  numOwnedCols++;
279  }
280 
281  // insert it into the new matrix
282  colIndices.resize(numOwnedCols);
283  colValues.resize(numOwnedCols);
284  mat->insertGlobalValues(contigRow, Teuchos::ArrayView<GO>(colIndices),
285  Teuchos::ArrayView<ST>(colValues));
286  colIndices.resize(maxNumEntries);
287  colValues.resize(maxNumEntries);
288  }
289 
290  // fill it and automagically optimize the storage
291  mat->fillComplete(colMap, rowMap);
292 
293  return mat;
294 }
295 
296 // build a single subblock Epetra_CrsMatrix
297 void rebuildSubBlock(int i, int j, const RCP<const Tpetra::CrsMatrix<ST, LO, GO, NT> >& A,
298  const std::vector<MapPair>& subMaps, Tpetra::CrsMatrix<ST, LO, GO, NT>& mat) {
299  // get the number of variables families
300  int numVarFamily = subMaps.size();
301 
302  TEUCHOS_ASSERT(i >= 0 && i < numVarFamily);
303  TEUCHOS_ASSERT(j >= 0 && j < numVarFamily);
304 
305  const Tpetra::Map<LO, GO, NT>& gRowMap = *subMaps[i].first; // new GIDs
306  const Tpetra::Map<LO, GO, NT>& rowMap = *subMaps[i].second; // contiguous GIDs
307  const Tpetra::Map<LO, GO, NT>& colMap = *subMaps[j].second;
308 
309  const RCP<Tpetra::Vector<GO, LO, GO, NT> > plocal2ContigGIDs =
310  getSubBlockColumnGIDs(*A, subMaps[j]);
311  Tpetra::Vector<GO, LO, GO, NT>& local2ContigGIDs = *plocal2ContigGIDs;
312 
313  mat.resumeFill();
314  mat.setAllToScalar(0.0);
315 
316  // get entry information
317  LO numMyRows = rowMap.getLocalNumElements();
318  LO maxNumEntries = A->getLocalMaxNumRowEntries();
319 
320  // for extraction
321  auto indices = typename Tpetra::CrsMatrix<ST, LO, GO, NT>::nonconst_local_inds_host_view_type(
322  Kokkos::ViewAllocateWithoutInitializing("rowIndices"), maxNumEntries);
323  auto values = typename Tpetra::CrsMatrix<ST, LO, GO, NT>::nonconst_values_host_view_type(
324  Kokkos::ViewAllocateWithoutInitializing("rowIndices"), maxNumEntries);
325 
326  // for insertion
327  std::vector<GO> colIndices(maxNumEntries);
328  std::vector<ST> colValues(maxNumEntries);
329 
330  const size_t invalid = Teuchos::OrdinalTraits<size_t>::invalid();
331 
332  // insert each row into subblock
333  // let FillComplete handle column distribution
334  for (LO localRow = 0; localRow < numMyRows; localRow++) {
335  size_t numEntries = invalid;
336  GO globalRow = gRowMap.getGlobalElement(localRow);
337  LO lid = A->getRowMap()->getLocalElement(globalRow);
338  GO contigRow = rowMap.getGlobalElement(localRow);
339  TEUCHOS_ASSERT(lid > -1);
340 
341  A->getLocalRowCopy(lid, indices, values, numEntries);
342 
343  LO numOwnedCols = 0;
344  for (size_t localCol = 0; localCol < numEntries; localCol++) {
345  TEUCHOS_ASSERT(indices(localCol) > -1);
346 
347  // if global id is not owned by this column
348  GO gid = local2ContigGIDs.getData()[indices(localCol)];
349  if (gid == -1) continue; // in contiguous row
350 
351  colIndices[numOwnedCols] = gid;
352  colValues[numOwnedCols] = values(localCol);
353  numOwnedCols++;
354  }
355 
356  // insert it into the new matrix
357  colIndices.resize(numOwnedCols);
358  colValues.resize(numOwnedCols);
359  mat.sumIntoGlobalValues(contigRow, Teuchos::ArrayView<GO>(colIndices),
360  Teuchos::ArrayView<ST>(colValues));
361  colIndices.resize(maxNumEntries);
362  colValues.resize(maxNumEntries);
363  }
364  mat.fillComplete(rcpFromRef(colMap), rcpFromRef(rowMap));
365 }
366 
367 } // namespace Blocking
368 } // namespace TpetraHelpers
369 } // namespace Teko
Teko_Utilities.hpp
Generated on Fri Nov 22 2024 09:19:41 for Teko by   doxygen 1.8.5