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