Teko  Version of the Day
 All Classes Files Functions Variables Pages
Teko_BlockingTpetra.cpp
1 /*
2 // @HEADER
3 //
4 // ***********************************************************************
5 //
6 // Teko: A package for block and physics based preconditioning
7 // Copyright 2010 Sandia Corporation
8 //
9 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
10 // the U.S. Government retains certain rights in this software.
11 //
12 // Redistribution and use in source and binary forms, with or without
13 // modification, are permitted provided that the following conditions are
14 // met:
15 //
16 // 1. Redistributions of source code must retain the above copyright
17 // notice, this list of conditions and the following disclaimer.
18 //
19 // 2. Redistributions in binary form must reproduce the above copyright
20 // notice, this list of conditions and the following disclaimer in the
21 // documentation and/or other materials provided with the distribution.
22 //
23 // 3. Neither the name of the Corporation nor the names of the
24 // contributors may be used to endorse or promote products derived from
25 // this software without specific prior written permission.
26 //
27 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38 //
39 // Questions? Contact Eric C. Cyr (eccyr@sandia.gov)
40 //
41 // ***********************************************************************
42 //
43 // @HEADER
44 
45 */
46 
47 #include "Teko_BlockingTpetra.hpp"
48 #include "Teko_Utilities.hpp"
49 
50 #include "Tpetra_Vector.hpp"
51 #include "Tpetra_Map.hpp"
52 
53 using Teuchos::RCP;
54 using Teuchos::rcp;
55 
56 namespace Teko {
57 namespace TpetraHelpers {
58 namespace Blocking {
59 
73 const MapPair buildSubMap(const std::vector<GO>& gid, const Teuchos::Comm<int>& comm) {
74  using GST = Tpetra::global_size_t;
75  const GST invalid = Teuchos::OrdinalTraits<GST>::invalid();
76  Teuchos::RCP<Tpetra::Map<LO, GO, NT> > gidMap = rcp(
77  new Tpetra::Map<LO, GO, NT>(invalid, Teuchos::ArrayView<const GO>(gid), 0, rcpFromRef(comm)));
78  Teuchos::RCP<Tpetra::Map<LO, GO, NT> > contigMap =
79  rcp(new Tpetra::Map<LO, GO, NT>(invalid, gid.size(), 0, rcpFromRef(comm)));
80 
81  return std::make_pair(gidMap, contigMap);
82 }
83 
93 const ImExPair buildExportImport(const Tpetra::Map<LO, GO, NT>& baseMap, const MapPair& maps) {
94  return std::make_pair(rcp(new Tpetra::Import<LO, GO, NT>(rcpFromRef(baseMap), maps.first)),
95  rcp(new Tpetra::Export<LO, GO, NT>(maps.first, rcpFromRef(baseMap))));
96 }
97 
105 void buildSubVectors(const std::vector<MapPair>& maps,
106  std::vector<RCP<Tpetra::MultiVector<ST, LO, GO, NT> > >& vectors, int count) {
107  std::vector<MapPair>::const_iterator mapItr;
108 
109  // loop over all maps
110  for (mapItr = maps.begin(); mapItr != maps.end(); mapItr++) {
111  // add new elements to vectors
112  RCP<Tpetra::MultiVector<ST, LO, GO, NT> > mv =
113  rcp(new Tpetra::MultiVector<ST, LO, GO, NT>((*mapItr).second, count));
114  vectors.push_back(mv);
115  }
116 }
117 
125 void one2many(std::vector<RCP<Tpetra::MultiVector<ST, LO, GO, NT> > >& many,
126  const Tpetra::MultiVector<ST, LO, GO, NT>& single,
127  const std::vector<RCP<Tpetra::Import<LO, GO, NT> > >& subImport) {
128  // std::vector<RCP<Epetra_Vector> >::const_iterator vecItr;
129  std::vector<RCP<Tpetra::MultiVector<ST, LO, GO, NT> > >::const_iterator vecItr;
130  std::vector<RCP<Tpetra::Import<LO, GO, NT> > >::const_iterator impItr;
131 
132  // using Importers fill the sub vectors from the mama vector
133  for (vecItr = many.begin(), impItr = subImport.begin(); vecItr != many.end();
134  ++vecItr, ++impItr) {
135  // for ease of access to the destination
136  RCP<Tpetra::MultiVector<ST, LO, GO, NT> > destVec = *vecItr;
137 
138  // extract the map with global indicies from the current vector
139  const Tpetra::Map<LO, GO, NT>& globalMap = *(*impItr)->getTargetMap();
140 
141  // build the import vector as a view on the destination
142  GO lda = destVec->getStride();
143  GO destSize = destVec->getGlobalLength() * destVec->getNumVectors();
144  std::vector<ST> destArray(destSize);
145  Teuchos::ArrayView<ST> destVals(destArray);
146  destVec->get1dCopy(destVals, lda);
147  Tpetra::MultiVector<ST, LO, GO, NT> importVector(rcpFromRef(globalMap), destVals, lda,
148  destVec->getNumVectors());
149 
150  // perform the import
151  importVector.doImport(single, **impItr, Tpetra::INSERT);
152 
153  Tpetra::Import<LO, GO, NT> importer(destVec->getMap(), destVec->getMap());
154  importVector.replaceMap(destVec->getMap());
155  destVec->doExport(importVector, importer, Tpetra::INSERT);
156  }
157 }
158 
169 void many2one(Tpetra::MultiVector<ST, LO, GO, NT>& one,
170  const std::vector<RCP<const Tpetra::MultiVector<ST, LO, GO, NT> > >& many,
171  const std::vector<RCP<Tpetra::Export<LO, GO, NT> > >& subExport) {
172  // std::vector<RCP<const Epetra_Vector> >::const_iterator vecItr;
173  std::vector<RCP<const Tpetra::MultiVector<ST, LO, GO, NT> > >::const_iterator vecItr;
174  std::vector<RCP<Tpetra::Export<LO, GO, NT> > >::const_iterator expItr;
175 
176  // using Exporters fill the empty vector from the sub-vectors
177  for (vecItr = many.begin(), expItr = subExport.begin(); vecItr != many.end();
178  ++vecItr, ++expItr) {
179  // for ease of access to the source
180  RCP<const Tpetra::MultiVector<ST, LO, GO, NT> > srcVec = *vecItr;
181 
182  // extract the map with global indicies from the current vector
183  const Tpetra::Map<LO, GO, NT>& globalMap = *(*expItr)->getSourceMap();
184 
185  // build the export vector as a view of the destination
186  GO lda = srcVec->getStride();
187  GO srcSize = srcVec->getGlobalLength() * srcVec->getNumVectors();
188  std::vector<ST> srcArray(srcSize);
189  Teuchos::ArrayView<ST> srcVals(srcArray);
190  srcVec->get1dCopy(srcVals, lda);
191  Tpetra::MultiVector<ST, LO, GO, NT> exportVector(rcpFromRef(globalMap), srcVals, lda,
192  srcVec->getNumVectors());
193 
194  // perform the export
195  one.doExport(exportVector, **expItr, Tpetra::INSERT);
196  }
197 }
198 
204 RCP<Tpetra::Vector<GO, LO, GO, NT> > getSubBlockColumnGIDs(
205  const Tpetra::CrsMatrix<ST, LO, GO, NT>& A, const MapPair& mapPair) {
206  RCP<const Tpetra::Map<LO, GO, NT> > blkGIDMap = mapPair.first;
207  RCP<const Tpetra::Map<LO, GO, NT> > blkContigMap = mapPair.second;
208 
209  // fill index vector for rows
210  Tpetra::Vector<GO, LO, GO, NT> rIndex(A.getRowMap(), true);
211  for (size_t i = 0; i < A.getLocalNumRows(); i++) {
212  // LID is need to get to contiguous map
213  LO lid = blkGIDMap->getLocalElement(
214  A.getRowMap()->getGlobalElement(i)); // this LID makes me nervous
215  if (lid > -1)
216  rIndex.replaceLocalValue(i, blkContigMap->getGlobalElement(lid));
217  else
218  rIndex.replaceLocalValue(i, -1);
219  }
220 
221  // get relavant column indices
222 
223  Tpetra::Import<LO, GO, NT> import(A.getRowMap(), A.getColMap());
224  RCP<Tpetra::Vector<GO, LO, GO, NT> > cIndex =
225  rcp(new Tpetra::Vector<GO, LO, GO, NT>(A.getColMap(), true));
226  cIndex->doImport(rIndex, import, Tpetra::INSERT);
227 
228  return cIndex;
229 }
230 
231 // build a single subblock Epetra_CrsMatrix
232 RCP<Tpetra::CrsMatrix<ST, LO, GO, NT> > buildSubBlock(
233  int i, int j, const RCP<const Tpetra::CrsMatrix<ST, LO, GO, NT> >& A,
234  const std::vector<MapPair>& subMaps) {
235  // get the number of variables families
236  int numVarFamily = subMaps.size();
237 
238  TEUCHOS_ASSERT(i >= 0 && i < numVarFamily);
239  TEUCHOS_ASSERT(j >= 0 && j < numVarFamily);
240  const RCP<const Tpetra::Map<LO, GO, NT> > gRowMap = subMaps[i].first; // new GIDs
241  const RCP<const Tpetra::Map<LO, GO, NT> > rowMap = subMaps[i].second; // contiguous GIDs
242  const RCP<const Tpetra::Map<LO, GO, NT> > colMap = subMaps[j].second;
243 
244  const RCP<Tpetra::Vector<GO, LO, GO, NT> > plocal2ContigGIDs =
245  getSubBlockColumnGIDs(*A, subMaps[j]);
246  Tpetra::Vector<GO, LO, GO, NT>& local2ContigGIDs = *plocal2ContigGIDs;
247 
248  // get entry information
249  LO numMyRows = rowMap->getLocalNumElements();
250  LO maxNumEntries = A->getLocalMaxNumRowEntries();
251 
252  // for extraction
253  auto indices = typename Tpetra::CrsMatrix<ST, LO, GO, NT>::nonconst_local_inds_host_view_type(
254  Kokkos::ViewAllocateWithoutInitializing("rowIndices"), maxNumEntries);
255  auto values = typename Tpetra::CrsMatrix<ST, LO, GO, NT>::nonconst_values_host_view_type(
256  Kokkos::ViewAllocateWithoutInitializing("rowIndices"), maxNumEntries);
257 
258  // for insertion
259  std::vector<GO> colIndices(maxNumEntries);
260  std::vector<ST> colValues(maxNumEntries);
261 
262  // for counting
263  std::vector<size_t> nEntriesPerRow(numMyRows);
264 
265  const size_t invalid = Teuchos::OrdinalTraits<size_t>::invalid();
266 
267  // insert each row into subblock
268  // Count the number of entries per row in the new matrix
269  for (LO localRow = 0; localRow < numMyRows; localRow++) {
270  size_t numEntries = invalid;
271  GO globalRow = gRowMap->getGlobalElement(localRow);
272  LO lid = A->getRowMap()->getLocalElement(globalRow);
273  TEUCHOS_ASSERT(lid > -1);
274 
275  A->getLocalRowCopy(lid, indices, values, numEntries);
276 
277  LO numOwnedCols = 0;
278  for (size_t localCol = 0; localCol < numEntries; localCol++) {
279  TEUCHOS_ASSERT(indices(localCol) > -1);
280 
281  // if global id is not owned by this column
282 
283  GO gid = local2ContigGIDs.getData()[indices[localCol]];
284  if (gid == -1) continue; // in contiguous row
285  numOwnedCols++;
286  }
287  nEntriesPerRow[localRow] = numOwnedCols;
288  }
289 
290  RCP<Tpetra::CrsMatrix<ST, LO, GO, NT> > mat = rcp(
291  new Tpetra::CrsMatrix<ST, LO, GO, NT>(rowMap, Teuchos::ArrayView<size_t>(nEntriesPerRow)));
292 
293  // insert each row into subblock
294  // let FillComplete handle column distribution
295  for (LO localRow = 0; localRow < numMyRows; localRow++) {
296  size_t numEntries = invalid;
297  GO globalRow = gRowMap->getGlobalElement(localRow);
298  LO lid = A->getRowMap()->getLocalElement(globalRow);
299  GO contigRow = rowMap->getGlobalElement(localRow);
300  TEUCHOS_ASSERT(lid > -1);
301 
302  A->getLocalRowCopy(lid, indices, values, numEntries);
303 
304  LO numOwnedCols = 0;
305  for (size_t localCol = 0; localCol < numEntries; localCol++) {
306  TEUCHOS_ASSERT(indices(localCol) > -1);
307 
308  // if global id is not owned by this column
309 
310  GO gid = local2ContigGIDs.getData()[indices(localCol)];
311  if (gid == -1) continue; // in contiguous row
312 
313  colIndices[numOwnedCols] = gid;
314  colValues[numOwnedCols] = values(localCol);
315  numOwnedCols++;
316  }
317 
318  // insert it into the new matrix
319  colIndices.resize(numOwnedCols);
320  colValues.resize(numOwnedCols);
321  mat->insertGlobalValues(contigRow, Teuchos::ArrayView<GO>(colIndices),
322  Teuchos::ArrayView<ST>(colValues));
323  colIndices.resize(maxNumEntries);
324  colValues.resize(maxNumEntries);
325  }
326 
327  // fill it and automagically optimize the storage
328  mat->fillComplete(colMap, rowMap);
329 
330  return mat;
331 }
332 
333 // build a single subblock Epetra_CrsMatrix
334 void rebuildSubBlock(int i, int j, const RCP<const Tpetra::CrsMatrix<ST, LO, GO, NT> >& A,
335  const std::vector<MapPair>& subMaps, Tpetra::CrsMatrix<ST, LO, GO, NT>& mat) {
336  // get the number of variables families
337  int numVarFamily = subMaps.size();
338 
339  TEUCHOS_ASSERT(i >= 0 && i < numVarFamily);
340  TEUCHOS_ASSERT(j >= 0 && j < numVarFamily);
341 
342  const Tpetra::Map<LO, GO, NT>& gRowMap = *subMaps[i].first; // new GIDs
343  const Tpetra::Map<LO, GO, NT>& rowMap = *subMaps[i].second; // contiguous GIDs
344  const Tpetra::Map<LO, GO, NT>& colMap = *subMaps[j].second;
345 
346  const RCP<Tpetra::Vector<GO, LO, GO, NT> > plocal2ContigGIDs =
347  getSubBlockColumnGIDs(*A, subMaps[j]);
348  Tpetra::Vector<GO, LO, GO, NT>& local2ContigGIDs = *plocal2ContigGIDs;
349 
350  mat.resumeFill();
351  mat.setAllToScalar(0.0);
352 
353  // get entry information
354  LO numMyRows = rowMap.getLocalNumElements();
355  LO maxNumEntries = A->getLocalMaxNumRowEntries();
356 
357  // for extraction
358  auto indices = typename Tpetra::CrsMatrix<ST, LO, GO, NT>::nonconst_local_inds_host_view_type(
359  Kokkos::ViewAllocateWithoutInitializing("rowIndices"), maxNumEntries);
360  auto values = typename Tpetra::CrsMatrix<ST, LO, GO, NT>::nonconst_values_host_view_type(
361  Kokkos::ViewAllocateWithoutInitializing("rowIndices"), maxNumEntries);
362 
363  // for insertion
364  std::vector<GO> colIndices(maxNumEntries);
365  std::vector<ST> colValues(maxNumEntries);
366 
367  const size_t invalid = Teuchos::OrdinalTraits<size_t>::invalid();
368 
369  // insert each row into subblock
370  // let FillComplete handle column distribution
371  for (LO localRow = 0; localRow < numMyRows; localRow++) {
372  size_t numEntries = invalid;
373  GO globalRow = gRowMap.getGlobalElement(localRow);
374  LO lid = A->getRowMap()->getLocalElement(globalRow);
375  GO contigRow = rowMap.getGlobalElement(localRow);
376  TEUCHOS_ASSERT(lid > -1);
377 
378  A->getLocalRowCopy(lid, indices, values, numEntries);
379 
380  LO numOwnedCols = 0;
381  for (size_t localCol = 0; localCol < numEntries; localCol++) {
382  TEUCHOS_ASSERT(indices(localCol) > -1);
383 
384  // if global id is not owned by this column
385  GO gid = local2ContigGIDs.getData()[indices(localCol)];
386  if (gid == -1) continue; // in contiguous row
387 
388  colIndices[numOwnedCols] = gid;
389  colValues[numOwnedCols] = values(localCol);
390  numOwnedCols++;
391  }
392 
393  // insert it into the new matrix
394  colIndices.resize(numOwnedCols);
395  colValues.resize(numOwnedCols);
396  mat.sumIntoGlobalValues(contigRow, Teuchos::ArrayView<GO>(colIndices),
397  Teuchos::ArrayView<ST>(colValues));
398  colIndices.resize(maxNumEntries);
399  colValues.resize(maxNumEntries);
400  }
401  mat.fillComplete(rcpFromRef(colMap), rcpFromRef(rowMap));
402 }
403 
404 } // namespace Blocking
405 } // namespace TpetraHelpers
406 } // namespace Teko
Teko_Utilities.hpp
Generated on Thu Apr 25 2024 09:20:43 for Teko by   doxygen 1.8.5