Teko Version of the Day
Loading...
Searching...
No Matches
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
53using Teuchos::RCP;
54using Teuchos::rcp;
55
56namespace Teko {
57namespace TpetraHelpers {
58namespace Blocking {
59
73const 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
92const 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
105void 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
124void 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
166void 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
201RCP<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.getLocalNumRows();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
227RCP<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->getLocalNumElements();
244 LO maxNumEntries = A->getLocalMaxNumRowEntries();
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
324void 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.getLocalNumElements();
344 LO maxNumEntries = A->getLocalMaxNumRowEntries();
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