Teko  Version of the Day
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 {
75  using GST = Tpetra::global_size_t;
76  const GST invalid = Teuchos::OrdinalTraits<GST>::invalid();
77  Teuchos::RCP<Tpetra::Map<LO,GO,NT> > gidMap = rcp(new Tpetra::Map<LO,GO,NT>(invalid,Teuchos::ArrayView<const GO>(gid),0,rcpFromRef(comm)));
78  Teuchos::RCP<Tpetra::Map<LO,GO,NT> > contigMap = rcp(new Tpetra::Map<LO,GO,NT>(invalid,gid.size(),0,rcpFromRef(comm)));
79 
80  return std::make_pair(gidMap,contigMap);
81 }
82 
92 const ImExPair buildExportImport(const Tpetra::Map<LO,GO,NT> & baseMap,const MapPair & maps)
93 {
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 {
108  std::vector<MapPair>::const_iterator mapItr;
109 
110  // loop over all maps
111  for(mapItr=maps.begin();mapItr!=maps.end();mapItr++) {
112  // add new elements to vectors
113  RCP<Tpetra::MultiVector<ST,LO,GO,NT> > mv = rcp(new Tpetra::MultiVector<ST,LO,GO,NT>((*mapItr).second,count));
114  vectors.push_back(mv);
115  }
116 }
117 
124 void one2many(std::vector<RCP<Tpetra::MultiVector<ST,LO,GO,NT> > > & many, const Tpetra::MultiVector<ST,LO,GO,NT> & single,
125  const std::vector<RCP<Tpetra::Import<LO,GO,NT> > > & subImport)
126 {
127  // std::vector<RCP<Epetra_Vector> >::const_iterator vecItr;
128  std::vector<RCP<Tpetra::MultiVector<ST,LO,GO,NT> > >::const_iterator vecItr;
129  std::vector<RCP<Tpetra::Import<LO,GO,NT> > >::const_iterator impItr;
130 
131  // using Importers fill the sub vectors from the mama vector
132  for(vecItr=many.begin(),impItr=subImport.begin();
133  vecItr!=many.end();++vecItr,++impItr) {
134  // for ease of access to the destination
135  RCP<Tpetra::MultiVector<ST,LO,GO,NT> > destVec = *vecItr;
136 
137  // extract the map with global indicies from the current vector
138  const Tpetra::Map<LO,GO,NT> & globalMap = *(*impItr)->getTargetMap();
139 
140  // build the import vector as a view on the destination
141  GO lda = destVec->getStride();
142  GO destSize = destVec->getGlobalLength()*destVec->getNumVectors();
143  std::vector<ST> destArray(destSize);
144  Teuchos::ArrayView<ST> destVals(destArray);
145  destVec->get1dCopy(destVals,lda);
146  Tpetra::MultiVector<ST,LO,GO,NT> importVector(rcpFromRef(globalMap),destVals,lda,destVec->getNumVectors());
147 
148  // perform the import
149  importVector.doImport(single,**impItr,Tpetra::INSERT);
150 
151  Tpetra::Import<LO,GO,NT> importer(destVec->getMap(),destVec->getMap());
152  importVector.replaceMap(destVec->getMap());
153  destVec->doExport(importVector,importer,Tpetra::INSERT);
154  }
155 }
156 
166 void many2one(Tpetra::MultiVector<ST,LO,GO,NT> & one, const std::vector<RCP<const Tpetra::MultiVector<ST,LO,GO,NT> > > & many,
167  const std::vector<RCP<Tpetra::Export<LO,GO,NT> > > & subExport)
168 {
169  // std::vector<RCP<const Epetra_Vector> >::const_iterator vecItr;
170  std::vector<RCP<const Tpetra::MultiVector<ST,LO,GO,NT> > >::const_iterator vecItr;
171  std::vector<RCP<Tpetra::Export<LO,GO,NT> > >::const_iterator expItr;
172 
173  // using Exporters fill the empty vector from the sub-vectors
174  for(vecItr=many.begin(),expItr=subExport.begin();
175  vecItr!=many.end();++vecItr,++expItr) {
176 
177  // for ease of access to the source
178  RCP<const Tpetra::MultiVector<ST,LO,GO,NT> > srcVec = *vecItr;
179 
180  // extract the map with global indicies from the current vector
181  const Tpetra::Map<LO,GO,NT> & globalMap = *(*expItr)->getSourceMap();
182 
183  // build the export vector as a view of the destination
184  GO lda = srcVec->getStride();
185  GO srcSize = srcVec->getGlobalLength()*srcVec->getNumVectors();
186  std::vector<ST> srcArray(srcSize);
187  Teuchos::ArrayView<ST> srcVals(srcArray);
188  srcVec->get1dCopy(srcVals,lda);
189  Tpetra::MultiVector<ST,LO,GO,NT> exportVector(rcpFromRef(globalMap),srcVals,lda,srcVec->getNumVectors());
190 
191  // perform the export
192  one.doExport(exportVector,**expItr,Tpetra::INSERT);
193  }
194 }
195 
201 RCP<Tpetra::Vector<GO,LO,GO,NT> > getSubBlockColumnGIDs(const Tpetra::CrsMatrix<ST,LO,GO,NT> & A,const MapPair & mapPair)
202 {
203  RCP<const Tpetra::Map<LO,GO,NT> > blkGIDMap = mapPair.first;
204  RCP<const Tpetra::Map<LO,GO,NT> > blkContigMap = mapPair.second;
205 
206  // fill index vector for rows
207  Tpetra::Vector<GO,LO,GO,NT> rIndex(A.getRowMap(),true);
208  for(size_t i=0;i<A.getNodeNumRows();i++) {
209  // LID is need to get to contiguous map
210  LO lid = blkGIDMap->getLocalElement(A.getRowMap()->getGlobalElement(i)); // this LID makes me nervous
211  if(lid>-1)
212  rIndex.replaceLocalValue(i,blkContigMap->getGlobalElement(lid));
213  else
214  rIndex.replaceLocalValue(i,-1);
215  }
216 
217  // get relavant column indices
218 
219  Tpetra::Import<LO,GO,NT> import(A.getRowMap(),A.getColMap());
220  RCP<Tpetra::Vector<GO,LO,GO,NT> > cIndex = rcp(new Tpetra::Vector<GO,LO,GO,NT>(A.getColMap(),true));
221  cIndex->doImport(rIndex,import,Tpetra::INSERT);
222 
223  return cIndex;
224 }
225 
226 // build a single subblock Epetra_CrsMatrix
227 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > buildSubBlock(int i,int j,const RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> >& A,const std::vector<MapPair> & subMaps)
228 {
229  // get the number of variables families
230  int numVarFamily = subMaps.size();
231 
232  TEUCHOS_ASSERT(i>=0 && i<numVarFamily);
233  TEUCHOS_ASSERT(j>=0 && j<numVarFamily);
234  const RCP<const Tpetra::Map<LO,GO,NT> > gRowMap = subMaps[i].first; // new GIDs
235  const RCP<const Tpetra::Map<LO,GO,NT> > rowMap = subMaps[i].second; // contiguous GIDs
236  const RCP<const Tpetra::Map<LO,GO,NT> > colMap = subMaps[j].second;
237 
238  const RCP<Tpetra::Vector<GO,LO,GO,NT> > plocal2ContigGIDs = getSubBlockColumnGIDs(*A,subMaps[j]);
239  Tpetra::Vector<GO,LO,GO,NT> & local2ContigGIDs = *plocal2ContigGIDs;
240 
241 
242  // get entry information
243  LO numMyRows = rowMap->getNodeNumElements();
244  LO maxNumEntries = A->getNodeMaxNumRowEntries();
245 
246  // for extraction
247  auto indices = typename Tpetra::CrsMatrix<ST,LO,GO,NT>::nonconst_local_inds_host_view_type(Kokkos::ViewAllocateWithoutInitializing("rowIndices"),maxNumEntries);
248  auto values = typename Tpetra::CrsMatrix<ST,LO,GO,NT>::nonconst_values_host_view_type(Kokkos::ViewAllocateWithoutInitializing("rowIndices"),maxNumEntries);
249 
250  // for insertion
251  std::vector<GO> colIndices(maxNumEntries);
252  std::vector<ST> colValues(maxNumEntries);
253 
254  // for counting
255  std::vector<size_t> nEntriesPerRow(numMyRows);
256 
257  const size_t invalid = Teuchos::OrdinalTraits<size_t>::invalid();
258 
259  // insert each row into subblock
260  // Count the number of entries per row in the new matrix
261  for(LO localRow=0;localRow<numMyRows;localRow++) {
262  size_t numEntries = invalid;
263  GO globalRow = gRowMap->getGlobalElement(localRow);
264  LO lid = A->getRowMap()->getLocalElement(globalRow);
265  TEUCHOS_ASSERT(lid>-1);
266 
267  A->getLocalRowCopy(lid, indices, values, numEntries);
268 
269  LO numOwnedCols = 0;
270  for(size_t localCol=0;localCol<numEntries;localCol++) {
271  TEUCHOS_ASSERT(indices(localCol)>-1);
272 
273  // if global id is not owned by this column
274 
275  GO gid = local2ContigGIDs.getData()[indices[localCol]];
276  if(gid==-1) continue; // in contiguous row
277  numOwnedCols++;
278  }
279  nEntriesPerRow[localRow] = numOwnedCols;
280  }
281 
282  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > mat = rcp(new Tpetra::CrsMatrix<ST,LO,GO,NT>(rowMap,Teuchos::ArrayView<size_t>(nEntriesPerRow)));
283 
284  // insert each row into subblock
285  // let FillComplete handle column distribution
286  for(LO localRow=0;localRow<numMyRows;localRow++) {
287  size_t numEntries = invalid;
288  GO globalRow = gRowMap->getGlobalElement(localRow);
289  LO lid = A->getRowMap()->getLocalElement(globalRow);
290  GO contigRow = rowMap->getGlobalElement(localRow);
291  TEUCHOS_ASSERT(lid>-1);
292 
293  A->getLocalRowCopy(lid, indices, values, numEntries);
294 
295  LO numOwnedCols = 0;
296  for(size_t localCol=0;localCol<numEntries;localCol++) {
297  TEUCHOS_ASSERT(indices(localCol)>-1);
298 
299  // if global id is not owned by this column
300 
301  GO gid = local2ContigGIDs.getData()[indices(localCol)];
302  if(gid==-1) continue; // in contiguous row
303 
304  colIndices[numOwnedCols] = gid;
305  colValues[numOwnedCols] = values(localCol);
306  numOwnedCols++;
307  }
308 
309  // insert it into the new matrix
310  colIndices.resize(numOwnedCols);
311  colValues.resize(numOwnedCols);
312  mat->insertGlobalValues(contigRow,Teuchos::ArrayView<GO>(colIndices),Teuchos::ArrayView<ST>(colValues));
313  colIndices.resize(maxNumEntries);
314  colValues.resize(maxNumEntries);
315  }
316 
317  // fill it and automagically optimize the storage
318  mat->fillComplete(colMap,rowMap);
319 
320  return mat;
321 }
322 
323 // build a single subblock Epetra_CrsMatrix
324 void rebuildSubBlock(int i,int j,const RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> >& A,const std::vector<MapPair> & subMaps,Tpetra::CrsMatrix<ST,LO,GO,NT> & mat)
325 {
326  // get the number of variables families
327  int numVarFamily = subMaps.size();
328 
329  TEUCHOS_ASSERT(i>=0 && i<numVarFamily);
330  TEUCHOS_ASSERT(j>=0 && j<numVarFamily);
331 
332  const Tpetra::Map<LO,GO,NT> & gRowMap = *subMaps[i].first; // new GIDs
333  const Tpetra::Map<LO,GO,NT> & rowMap = *subMaps[i].second; // contiguous GIDs
334  const Tpetra::Map<LO,GO,NT> & colMap = *subMaps[j].second;
335 
336  const RCP<Tpetra::Vector<GO,LO,GO,NT> > plocal2ContigGIDs = getSubBlockColumnGIDs(*A,subMaps[j]);
337  Tpetra::Vector<GO,LO,GO,NT> & local2ContigGIDs = *plocal2ContigGIDs;
338 
339  mat.resumeFill();
340  mat.setAllToScalar(0.0);
341 
342  // get entry information
343  LO numMyRows = rowMap.getNodeNumElements();
344  LO maxNumEntries = A->getNodeMaxNumRowEntries();
345 
346  // for extraction
347  auto indices = typename Tpetra::CrsMatrix<ST,LO,GO,NT>::nonconst_local_inds_host_view_type(Kokkos::ViewAllocateWithoutInitializing("rowIndices"),maxNumEntries);
348  auto values = typename Tpetra::CrsMatrix<ST,LO,GO,NT>::nonconst_values_host_view_type(Kokkos::ViewAllocateWithoutInitializing("rowIndices"),maxNumEntries);
349 
350  // for insertion
351  std::vector<GO> colIndices(maxNumEntries);
352  std::vector<ST> colValues(maxNumEntries);
353 
354  const size_t invalid = Teuchos::OrdinalTraits<size_t>::invalid();
355 
356  // insert each row into subblock
357  // let FillComplete handle column distribution
358  for(LO localRow=0;localRow<numMyRows;localRow++) {
359  size_t numEntries = invalid;
360  GO globalRow = gRowMap.getGlobalElement(localRow);
361  LO lid = A->getRowMap()->getLocalElement(globalRow);
362  GO contigRow = rowMap.getGlobalElement(localRow);
363  TEUCHOS_ASSERT(lid>-1);
364 
365  A->getLocalRowCopy(lid, indices, values, numEntries);
366 
367  LO numOwnedCols = 0;
368  for(size_t localCol=0;localCol<numEntries;localCol++) {
369  TEUCHOS_ASSERT(indices(localCol)>-1);
370 
371  // if global id is not owned by this column
372  GO gid = local2ContigGIDs.getData()[indices(localCol)];
373  if(gid==-1) continue; // in contiguous row
374 
375  colIndices[numOwnedCols] = gid;
376  colValues[numOwnedCols] = values(localCol);
377  numOwnedCols++;
378  }
379 
380  // insert it into the new matrix
381  colIndices.resize(numOwnedCols);
382  colValues.resize(numOwnedCols);
383  mat.sumIntoGlobalValues(contigRow,Teuchos::ArrayView<GO>(colIndices),Teuchos::ArrayView<ST>(colValues));
384  colIndices.resize(maxNumEntries);
385  colValues.resize(maxNumEntries);
386  }
387  mat.fillComplete(rcpFromRef(colMap),rcpFromRef(rowMap));
388 }
389 
390 } // end Blocking
391 } // end Epetra
392 } // end Teko
Teko_Utilities.hpp
Generated by doxygen 1.9.1