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 
240  // get entry information
241  LO numMyRows = rowMap->getNodeNumElements();
242  LO maxNumEntries = A->getNodeMaxNumRowEntries();
243 
244  // for extraction
245  std::vector<LO> indices(maxNumEntries);
246  std::vector<ST> values(maxNumEntries);
247 
248  // for insertion
249  std::vector<GO> colIndices(maxNumEntries);
250  std::vector<ST> colValues(maxNumEntries);
251 
252  // for counting
253  std::vector<size_t> nEntriesPerRow(numMyRows);
254 
255  // insert each row into subblock
256  // Count the number of entries per row in the new matrix
257  for(LO localRow=0;localRow<numMyRows;localRow++) {
258  size_t numEntries = -1;
259  GO globalRow = gRowMap->getGlobalElement(localRow);
260  LO lid = A->getRowMap()->getLocalElement(globalRow);
261  TEUCHOS_ASSERT(lid>-1);
262 
263  A->getLocalRowCopy(lid, Teuchos::ArrayView<LO>(indices), Teuchos::ArrayView<ST>(values),numEntries);
264 
265  LO numOwnedCols = 0;
266  for(size_t localCol=0;localCol<numEntries;localCol++) {
267  TEUCHOS_ASSERT(indices[localCol]>-1);
268 
269  // if global id is not owned by this column
270 
271  GO gid = local2ContigGIDs.getData()[indices[localCol]];
272  if(gid==-1) continue; // in contiguous row
273  numOwnedCols++;
274  }
275  nEntriesPerRow[localRow] = numOwnedCols;
276  }
277 
278  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > mat = rcp(new Tpetra::CrsMatrix<ST,LO,GO,NT>(rowMap,Teuchos::ArrayView<size_t>(nEntriesPerRow)));
279 
280  // insert each row into subblock
281  // let FillComplete handle column distribution
282  for(LO localRow=0;localRow<numMyRows;localRow++) {
283  size_t numEntries = -1;
284  GO globalRow = gRowMap->getGlobalElement(localRow);
285  LO lid = A->getRowMap()->getLocalElement(globalRow);
286  GO contigRow = rowMap->getGlobalElement(localRow);
287  TEUCHOS_ASSERT(lid>-1);
288 
289  A->getLocalRowCopy(lid, Teuchos::ArrayView<LO>(indices), Teuchos::ArrayView<ST>(values),numEntries);
290 
291  LO numOwnedCols = 0;
292  for(size_t localCol=0;localCol<numEntries;localCol++) {
293  TEUCHOS_ASSERT(indices[localCol]>-1);
294 
295  // if global id is not owned by this column
296 
297  GO gid = local2ContigGIDs.getData()[indices[localCol]];
298  if(gid==-1) continue; // in contiguous row
299 
300  colIndices[numOwnedCols] = gid;
301  colValues[numOwnedCols] = values[localCol];
302  numOwnedCols++;
303  }
304 
305  // insert it into the new matrix
306  colIndices.resize(numOwnedCols);
307  colValues.resize(numOwnedCols);
308  mat->insertGlobalValues(contigRow,Teuchos::ArrayView<GO>(colIndices),Teuchos::ArrayView<ST>(colValues));
309  colIndices.resize(maxNumEntries);
310  colValues.resize(maxNumEntries);
311  }
312 
313  // fill it and automagically optimize the storage
314  mat->fillComplete(colMap,rowMap);
315 
316  return mat;
317 }
318 
319 // build a single subblock Epetra_CrsMatrix
320 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)
321 {
322  // get the number of variables families
323  int numVarFamily = subMaps.size();
324 
325  TEUCHOS_ASSERT(i>=0 && i<numVarFamily);
326  TEUCHOS_ASSERT(j>=0 && j<numVarFamily);
327 
328  const Tpetra::Map<LO,GO,NT> & gRowMap = *subMaps[i].first; // new GIDs
329  const Tpetra::Map<LO,GO,NT> & rowMap = *subMaps[i].second; // contiguous GIDs
330  const Tpetra::Map<LO,GO,NT> & colMap = *subMaps[j].second;
331 
332  const RCP<Tpetra::Vector<GO,LO,GO,NT> > plocal2ContigGIDs = getSubBlockColumnGIDs(*A,subMaps[j]);
333  Tpetra::Vector<GO,LO,GO,NT> & local2ContigGIDs = *plocal2ContigGIDs;
334 
335  mat.resumeFill();
336  mat.setAllToScalar(0.0);
337 
338  // get entry information
339  LO numMyRows = rowMap.getNodeNumElements();
340  LO maxNumEntries = A->getNodeMaxNumRowEntries();
341 
342  // for extraction
343  std::vector<LO> indices(maxNumEntries);
344  std::vector<ST> values(maxNumEntries);
345 
346  // for insertion
347  std::vector<GO> colIndices(maxNumEntries);
348  std::vector<ST> colValues(maxNumEntries);
349 
350  // insert each row into subblock
351  // let FillComplete handle column distribution
352  for(LO localRow=0;localRow<numMyRows;localRow++) {
353  size_t numEntries = -1;
354  GO globalRow = gRowMap.getGlobalElement(localRow);
355  LO lid = A->getRowMap()->getLocalElement(globalRow);
356  GO contigRow = rowMap.getGlobalElement(localRow);
357  TEUCHOS_ASSERT(lid>-1);
358 
359  A->getLocalRowCopy(lid, Teuchos::ArrayView<LO>(indices), Teuchos::ArrayView<ST>(values), numEntries);
360 
361  LO numOwnedCols = 0;
362  for(size_t localCol=0;localCol<numEntries;localCol++) {
363  TEUCHOS_ASSERT(indices[localCol]>-1);
364 
365  // if global id is not owned by this column
366  GO gid = local2ContigGIDs.getData()[indices[localCol]];
367  if(gid==-1) continue; // in contiguous row
368 
369  colIndices[numOwnedCols] = gid;
370  colValues[numOwnedCols] = values[localCol];
371  numOwnedCols++;
372  }
373 
374  // insert it into the new matrix
375  colIndices.resize(numOwnedCols);
376  colValues.resize(numOwnedCols);
377  mat.sumIntoGlobalValues(contigRow,Teuchos::ArrayView<GO>(colIndices),Teuchos::ArrayView<ST>(colValues));
378  colIndices.resize(maxNumEntries);
379  colValues.resize(maxNumEntries);
380  }
381  mat.fillComplete(rcpFromRef(colMap),rcpFromRef(rowMap));
382 }
383 
384 } // end Blocking
385 } // end Epetra
386 } // end Teko
Teko_Utilities.hpp
Generated by   doxygen 1.8.5