Teko  Version of the Day
Teko_InterlacedTpetra.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_InterlacedTpetra.hpp"
48 #include "Tpetra_Import.hpp"
49 
50 #include <vector>
51 
52 using Teuchos::RCP;
53 using Teuchos::rcp;
54 
55 namespace Teko {
56 namespace TpetraHelpers {
57 namespace Strided {
58 
59 // this assumes that there are numGlobals with numVars each interlaced
60 // i.e. for numVars = 2 (u,v) then the vector is
61 // [u_0,v_0,u_1,v_1,u_2,v_2, ..., u_(numGlobals-1),v_(numGlobals-1)]
62 void buildSubMaps(GO numGlobals,int numVars,const Teuchos::Comm<int> & comm,std::vector<std::pair<int,RCP<Tpetra::Map<LO,GO,NT> > > > & subMaps)
63 {
64  std::vector<int> vars;
65 
66  // build vector describing the sub maps
67  for(int i=0;i<numVars;i++) vars.push_back(1);
68 
69  // build all the submaps
70  buildSubMaps(numGlobals,vars,comm,subMaps);
71 }
72 
73 // build maps to make other conversions
74 void buildSubMaps(const Tpetra::Map<LO,GO,NT> & globalMap,const std::vector<int> & vars,const Teuchos::Comm<int> & comm,
75  std::vector<std::pair<int,Teuchos::RCP<Tpetra::Map<LO,GO,NT> > > > & subMaps)
76 {
77  buildSubMaps(globalMap.getGlobalNumElements(),globalMap.getNodeNumElements(),globalMap.getMinGlobalIndex(),
78  vars,comm,subMaps);
79 }
80 
81 // build maps to make other conversions
82 void buildSubMaps(GO numGlobals,const std::vector<int> & vars,const Teuchos::Comm<int> & comm,std::vector<std::pair<int,Teuchos::RCP<Tpetra::Map<LO,GO,NT> > > > & subMaps)
83 {
84  std::vector<int>::const_iterator varItr;
85 
86  // compute total number of variables
87  int numGlobalVars = 0;
88  for(varItr=vars.begin();varItr!=vars.end();++varItr)
89  numGlobalVars += *varItr;
90 
91  // must be an even number of globals
92  TEUCHOS_ASSERT((numGlobals%numGlobalVars)==0);
93 
94  Tpetra::Map<LO,GO,NT> sampleMap(numGlobals/numGlobalVars,0,rcpFromRef(comm));
95 
96  buildSubMaps(numGlobals,numGlobalVars*sampleMap.getNodeNumElements(),numGlobalVars*sampleMap.getMinGlobalIndex(),vars,comm,subMaps);
97 }
98 
99 // build maps to make other conversions
100 void buildSubMaps(GO numGlobals,LO numMyElements,GO minMyGID,const std::vector<int> & vars,const Teuchos::Comm<int> & comm,
101  std::vector<std::pair<int,Teuchos::RCP<Tpetra::Map<LO,GO,NT> > > > & subMaps)
102 {
103  std::vector<int>::const_iterator varItr;
104 
105  // compute total number of variables
106  int numGlobalVars = 0;
107  for(varItr=vars.begin();varItr!=vars.end();++varItr)
108  numGlobalVars += *varItr;
109 
110  // must be an even number of globals
111  TEUCHOS_ASSERT((numGlobals%numGlobalVars)==0);
112  TEUCHOS_ASSERT((numMyElements%numGlobalVars)==0);
113  TEUCHOS_ASSERT((minMyGID%numGlobalVars)==0);
114 
115  LO numBlocks = numMyElements/numGlobalVars;
116  GO minBlockID = minMyGID/numGlobalVars;
117 
118  subMaps.clear();
119 
120  // index into local block in strided map
121  GO blockOffset = 0;
122  for(varItr=vars.begin();varItr!=vars.end();++varItr) {
123  LO numLocalVars = *varItr;
124  GO numAllElmts = numLocalVars*numGlobals/numGlobalVars;
125 #ifndef NDEBUG
126  LO numMyElmts = numLocalVars * numBlocks;
127 #endif
128 
129  // create global arrays describing the as of yet uncreated maps
130  std::vector<GO> subGlobals;
131  std::vector<GO> contigGlobals; // the contiguous globals
132 
133  // loop over each block of variables
134  LO count = 0;
135  for(LO blockNum=0;blockNum<numBlocks;blockNum++) {
136 
137  // loop over each local variable in the block
138  for(LO local=0;local<numLocalVars;++local) {
139  // global block number = minGID+blockNum
140  // block begin global id = numGlobalVars*(minGID+blockNum)
141  // global id block offset = blockOffset+local
142  subGlobals.push_back((minBlockID+blockNum)*numGlobalVars+blockOffset+local);
143 
144  // also build the contiguous IDs
145  contigGlobals.push_back(numLocalVars*minBlockID+count);
146  count++;
147  }
148  }
149 
150  // sanity check
151  assert((size_t) numMyElmts==subGlobals.size());
152 
153  // create the map with contiguous elements and the map with global elements
154  RCP<Tpetra::Map<LO,GO,NT> > subMap = rcp(new Tpetra::Map<LO,GO,NT>(numAllElmts,Teuchos::ArrayView<GO>(subGlobals),0,rcpFromRef(comm)));
155  RCP<Tpetra::Map<LO,GO,NT> > contigMap = rcp(new Tpetra::Map<LO,GO,NT>(numAllElmts,Teuchos::ArrayView<GO>(contigGlobals),0,rcpFromRef(comm)));
156 
157  Teuchos::set_extra_data(contigMap,"contigMap",Teuchos::inOutArg(subMap));
158  subMaps.push_back(std::make_pair(numLocalVars,subMap));
159 
160  // update the block offset
161  blockOffset += numLocalVars;
162  }
163 }
164 
165 void buildExportImport(const Tpetra::Map<LO,GO,NT> & baseMap, const std::vector<std::pair<int,RCP<Tpetra::Map<LO,GO,NT> > > > & subMaps,
166  std::vector<RCP<Tpetra::Export<LO,GO,NT> > > & subExport,
167  std::vector<RCP<Tpetra::Import<LO,GO,NT> > > & subImport)
168 {
169  std::vector<std::pair<int,RCP<Tpetra::Map<LO,GO,NT> > > >::const_iterator mapItr;
170 
171  // build importers and exporters
172  for(mapItr=subMaps.begin();mapItr!=subMaps.end();++mapItr) {
173  // exctract basic map
174  const Tpetra::Map<LO,GO,NT> & map = *(mapItr->second);
175 
176  // add new elements to vectors
177  subImport.push_back(rcp(new Tpetra::Import<LO,GO,NT>(rcpFromRef(baseMap),rcpFromRef(map))));
178  subExport.push_back(rcp(new Tpetra::Export<LO,GO,NT>(rcpFromRef(map),rcpFromRef(baseMap))));
179  }
180 }
181 
182 void buildSubVectors(const std::vector<std::pair<int,RCP<Tpetra::Map<LO,GO,NT> > > > & subMaps,std::vector<RCP<Tpetra::MultiVector<ST,LO,GO,NT> > > & subVectors,int count)
183 {
184  std::vector<std::pair<int,RCP<Tpetra::Map<LO,GO,NT> > > >::const_iterator mapItr;
185 
186  // build vectors
187  for(mapItr=subMaps.begin();mapItr!=subMaps.end();++mapItr) {
188  // exctract basic map
189  const Tpetra::Map<LO,GO,NT> & map = *(Teuchos::get_extra_data<RCP<Tpetra::Map<LO,GO,NT> > >(mapItr->second,"contigMap"));
190 
191  // add new elements to vectors
192  RCP<Tpetra::MultiVector<ST,LO,GO,NT> > mv = rcp(new Tpetra::MultiVector<ST,LO,GO,NT>(rcpFromRef(map),count));
193  Teuchos::set_extra_data(mapItr->second,"globalMap",Teuchos::inOutArg(mv));
194  subVectors.push_back(mv);
195  }
196 }
197 
198 void associateSubVectors(const std::vector<std::pair<int,RCP<Tpetra::Map<LO,GO,NT> > > > & subMaps,std::vector<RCP<const Tpetra::MultiVector<ST,LO,GO,NT> > > & subVectors)
199 {
200  std::vector<std::pair<int,RCP<Tpetra::Map<LO,GO,NT> > > >::const_iterator mapItr;
201  std::vector<RCP<const Tpetra::MultiVector<ST,LO,GO,NT> > >::iterator vecItr;
202 
203  TEUCHOS_ASSERT(subMaps.size()==subVectors.size());
204 
205  // associate the sub vectors with the subMaps
206  for(mapItr=subMaps.begin(),vecItr=subVectors.begin();mapItr!=subMaps.end();++mapItr,++vecItr)
207  Teuchos::set_extra_data(mapItr->second,"globalMap",Teuchos::inOutArg(*vecItr),Teuchos::POST_DESTROY,false);
208 }
209 
210 // build a single subblock Epetra_CrsMatrix
211 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<std::pair<int,RCP<Tpetra::Map<LO,GO,NT> > > > & subMaps)
212 {
213  // get the number of variables families
214  int numVarFamily = subMaps.size();
215 
216  TEUCHOS_ASSERT(i>=0 && i<numVarFamily);
217  TEUCHOS_ASSERT(j>=0 && j<numVarFamily);
218 
219  const Tpetra::Map<LO,GO,NT> & gRowMap = *subMaps[i].second;
220  const Tpetra::Map<LO,GO,NT> & rowMap = *Teuchos::get_extra_data<RCP<Tpetra::Map<LO,GO,NT> > >(subMaps[i].second,"contigMap");
221  const Tpetra::Map<LO,GO,NT> & colMap = *Teuchos::get_extra_data<RCP<Tpetra::Map<LO,GO,NT> > >(subMaps[j].second,"contigMap");
222  int colFamilyCnt = subMaps[j].first;
223 
224  // compute the number of global variables
225  // and the row and column block offset
226  GO numGlobalVars = 0;
227  GO rowBlockOffset = 0;
228  GO colBlockOffset = 0;
229  for(int k=0;k<numVarFamily;k++) {
230  numGlobalVars += subMaps[k].first;
231 
232  // compute block offsets
233  if(k<i) rowBlockOffset += subMaps[k].first;
234  if(k<j) colBlockOffset += subMaps[k].first;
235  }
236 
237  // copy all global rows to here
238  Tpetra::Import<LO,GO,NT> import(A->getRowMap(),rcpFromRef(gRowMap));
239  Tpetra::CrsMatrix<ST,LO,GO,NT> localA(rcpFromRef(gRowMap),0);
240  localA.doImport(*A,import,Tpetra::INSERT);
241 
242  // get entry information
243  LO numMyRows = rowMap.getNodeNumElements();
244  LO maxNumEntries = A->getGlobalMaxNumRowEntries();
245 
246  // for extraction
247  auto indices = typename Tpetra::CrsMatrix<ST,LO,GO,NT>::nonconst_global_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 counting row sizes
251  std::vector<size_t> numEntriesPerRow(numMyRows,0);
252 
253  const size_t invalid = Teuchos::OrdinalTraits<size_t>::invalid();
254 
255  // Count the sizes of each row, using same logic as insertion below
256  for(LO localRow=0;localRow<numMyRows;localRow++) {
257  size_t numEntries = invalid;
258  GO globalRow = gRowMap.getGlobalElement(localRow);
259  GO contigRow = rowMap.getGlobalElement(localRow);
260 
261  TEUCHOS_ASSERT(globalRow>=0);
262  TEUCHOS_ASSERT(contigRow>=0);
263 
264  // extract a global row copy
265  localA.getGlobalRowCopy(globalRow, indices, values, numEntries);
266  LO numOwnedCols = 0;
267  for(size_t localCol=0;localCol<numEntries;localCol++) {
268  GO globalCol = indices(localCol);
269 
270  // determinate which block this column ID is in
271  int block = globalCol / numGlobalVars;
272 
273  bool inFamily = true;
274 
275  // test the beginning of the block
276  inFamily &= (block*numGlobalVars+colBlockOffset <= globalCol);
277  inFamily &= ((block*numGlobalVars+colBlockOffset+colFamilyCnt) > globalCol);
278 
279  // is this column in the variable family
280  if(inFamily) {
281  numOwnedCols++;
282  }
283  }
284  numEntriesPerRow[localRow] += numOwnedCols;
285  }
286 
287  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > mat =
288  rcp(new Tpetra::CrsMatrix<ST,LO,GO,NT>(rcpFromRef(rowMap), Teuchos::ArrayView<const size_t>(numEntriesPerRow)));
289 
290  // for insertion
291  std::vector<GO> colIndices(maxNumEntries);
292  std::vector<ST> colValues(maxNumEntries);
293 
294  // insert each row into subblock
295  // let FillComplete handle column distribution
296  for(LO localRow=0;localRow<numMyRows;localRow++) {
297  size_t numEntries = invalid;
298  GO globalRow = gRowMap.getGlobalElement(localRow);
299  GO contigRow = rowMap.getGlobalElement(localRow);
300 
301  TEUCHOS_ASSERT(globalRow>=0);
302  TEUCHOS_ASSERT(contigRow>=0);
303 
304  // extract a global row copy
305  localA.getGlobalRowCopy(globalRow, indices, values, numEntries);
306  LO numOwnedCols = 0;
307  for(size_t localCol=0;localCol<numEntries;localCol++) {
308  GO globalCol = indices(localCol);
309 
310  // determinate which block this column ID is in
311  int block = globalCol / numGlobalVars;
312 
313  bool inFamily = true;
314 
315  // test the beginning of the block
316  inFamily &= (block*numGlobalVars+colBlockOffset <= globalCol);
317  inFamily &= ((block*numGlobalVars+colBlockOffset+colFamilyCnt) > globalCol);
318 
319  // is this column in the variable family
320  if(inFamily) {
321  GO familyOffset = globalCol-(block*numGlobalVars+colBlockOffset);
322 
323  colIndices[numOwnedCols] = block*colFamilyCnt + familyOffset;
324  colValues[numOwnedCols] = values(localCol);
325 
326  numOwnedCols++;
327  }
328  }
329 
330  // insert it into the new matrix
331  colIndices.resize(numOwnedCols);
332  colValues.resize(numOwnedCols);
333  mat->insertGlobalValues(contigRow,Teuchos::ArrayView<GO>(colIndices),Teuchos::ArrayView<ST>(colValues));
334  colIndices.resize(maxNumEntries);
335  colValues.resize(maxNumEntries);
336  }
337 
338  // fill it and automagically optimize the storage
339  mat->fillComplete(rcpFromRef(colMap),rcpFromRef(rowMap));
340 
341  return mat;
342 }
343 
344 // rebuild a single subblock Epetra_CrsMatrix
345 void rebuildSubBlock(int i,int j,const RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > & A,const std::vector<std::pair<int,RCP<Tpetra::Map<LO,GO,NT> > > > & subMaps,Tpetra::CrsMatrix<ST,LO,GO,NT> & mat)
346 {
347  // get the number of variables families
348  int numVarFamily = subMaps.size();
349 
350  TEUCHOS_ASSERT(i>=0 && i<numVarFamily);
351  TEUCHOS_ASSERT(j>=0 && j<numVarFamily);
352  TEUCHOS_ASSERT(mat.isFillComplete());
353 
354  const Tpetra::Map<LO,GO,NT> & gRowMap = *subMaps[i].second;
355  const Tpetra::Map<LO,GO,NT> & rowMap = *Teuchos::get_extra_data<RCP<Tpetra::Map<LO,GO,NT> > >(subMaps[i].second,"contigMap");
356  const Tpetra::Map<LO,GO,NT> & colMap = *Teuchos::get_extra_data<RCP<Tpetra::Map<LO,GO,NT> > >(subMaps[j].second,"contigMap");
357  int colFamilyCnt = subMaps[j].first;
358 
359  // compute the number of global variables
360  // and the row and column block offset
361  GO numGlobalVars = 0;
362  GO rowBlockOffset = 0;
363  GO colBlockOffset = 0;
364  for(int k=0;k<numVarFamily;k++) {
365  numGlobalVars += subMaps[k].first;
366 
367  // compute block offsets
368  if(k<i) rowBlockOffset += subMaps[k].first;
369  if(k<j) colBlockOffset += subMaps[k].first;
370  }
371 
372  // copy all global rows to here
373  Tpetra::Import<LO,GO,NT> import(A->getRowMap(),rcpFromRef(gRowMap));
374  Tpetra::CrsMatrix<ST,LO,GO,NT> localA(rcpFromRef(gRowMap),0);
375  localA.doImport(*A,import,Tpetra::INSERT);
376 
377  // clear out the old matrix
378  mat.resumeFill();
379  mat.setAllToScalar(0.0);
380 
381  // get entry information
382  LO numMyRows = rowMap.getNodeNumElements();
383  GO maxNumEntries = A->getGlobalMaxNumRowEntries();
384 
385  // for extraction
386  auto indices = typename Tpetra::CrsMatrix<ST,LO,GO,NT>::nonconst_global_inds_host_view_type(Kokkos::ViewAllocateWithoutInitializing("rowIndices"),maxNumEntries);
387  auto values = typename Tpetra::CrsMatrix<ST,LO,GO,NT>::nonconst_values_host_view_type(Kokkos::ViewAllocateWithoutInitializing("rowIndices"),maxNumEntries);
388 
389  // for insertion
390  std::vector<GO> colIndices(maxNumEntries);
391  std::vector<ST> colValues(maxNumEntries);
392 
393  const size_t invalid = Teuchos::OrdinalTraits<size_t>::invalid();
394 
395  // insert each row into subblock
396  // let FillComplete handle column distribution
397  for(LO localRow=0;localRow<numMyRows;localRow++) {
398  size_t numEntries = invalid;
399  GO globalRow = gRowMap.getGlobalElement(localRow);
400  GO contigRow = rowMap.getGlobalElement(localRow);
401 
402  TEUCHOS_ASSERT(globalRow>=0);
403  TEUCHOS_ASSERT(contigRow>=0);
404 
405  // extract a global row copy
406  localA.getGlobalRowCopy(globalRow, indices, values, numEntries);
407 
408  LO numOwnedCols = 0;
409  for(size_t localCol=0;localCol<numEntries;localCol++) {
410  GO globalCol = indices(localCol);
411 
412  // determinate which block this column ID is in
413  int block = globalCol / numGlobalVars;
414 
415  bool inFamily = true;
416 
417  // test the beginning of the block
418  inFamily &= (block*numGlobalVars+colBlockOffset <= globalCol);
419  inFamily &= ((block*numGlobalVars+colBlockOffset+colFamilyCnt) > globalCol);
420 
421  // is this column in the variable family
422  if(inFamily) {
423  GO familyOffset = globalCol-(block*numGlobalVars+colBlockOffset);
424 
425  colIndices[numOwnedCols] = block*colFamilyCnt + familyOffset;
426  colValues[numOwnedCols] = values(localCol);
427 
428  numOwnedCols++;
429  }
430  }
431 
432  // insert it into the new matrix
433  colIndices.resize(numOwnedCols);
434  colValues.resize(numOwnedCols);
435  mat.sumIntoGlobalValues(contigRow,Teuchos::ArrayView<GO>(colIndices),Teuchos::ArrayView<ST>(colValues));
436  colIndices.resize(maxNumEntries);
437  colValues.resize(maxNumEntries);
438  }
439  mat.fillComplete(rcpFromRef(colMap),rcpFromRef(rowMap));
440 }
441 
442 
443 // collect subvectors into a single global vector
444 void many2one(Tpetra::MultiVector<ST,LO,GO,NT> & one, const std::vector<RCP<const Tpetra::MultiVector<ST,LO,GO,NT> > > & many,
445  const std::vector<RCP<Tpetra::Export<LO,GO,NT> > > & subExport)
446 {
447  // std::vector<RCP<const Epetra_Vector> >::const_iterator vecItr;
448  std::vector<RCP<const Tpetra::MultiVector<ST,LO,GO,NT> > >::const_iterator vecItr;
449  std::vector<RCP<Tpetra::Export<LO,GO,NT> > >::const_iterator expItr;
450 
451  // using Exporters fill the empty vector from the sub-vectors
452  for(vecItr=many.begin(),expItr=subExport.begin();
453  vecItr!=many.end();++vecItr,++expItr) {
454 
455  // for ease of access to the source
456  RCP<const Tpetra::MultiVector<ST,LO,GO,NT> > srcVec = *vecItr;
457 
458  // extract the map with global indicies from the current vector
459  const Tpetra::Map<LO,GO,NT> & globalMap = *(Teuchos::get_extra_data<RCP<Tpetra::Map<LO,GO,NT> > >(srcVec,"globalMap"));
460 
461  // build the export vector as a view of the destination
462  GO lda = srcVec->getStride();
463  GO srcSize = srcVec->getGlobalLength()*srcVec->getNumVectors();
464  std::vector<ST> srcArray(srcSize);
465  Teuchos::ArrayView<ST> srcVals(srcArray);
466  srcVec->get1dCopy(srcVals,lda);
467  Tpetra::MultiVector<ST,LO,GO,NT> exportVector(rcpFromRef(globalMap),srcVals,lda,srcVec->getNumVectors());
468 
469  // perform the export
470  one.doExport(exportVector,**expItr,Tpetra::INSERT);
471  }
472 }
473 
474 // distribute one global vector into a many subvectors
475 void one2many(std::vector<RCP<Tpetra::MultiVector<ST,LO,GO,NT> > > & many,const Tpetra::MultiVector<ST,LO,GO,NT> & single,
476  const std::vector<RCP<Tpetra::Import<LO,GO,NT> > > & subImport)
477 {
478  // std::vector<RCP<Epetra_Vector> >::const_iterator vecItr;
479  std::vector<RCP<Tpetra::MultiVector<ST,LO,GO,NT> > >::const_iterator vecItr;
480  std::vector<RCP<Tpetra::Import<LO,GO,NT> > >::const_iterator impItr;
481 
482  // using Importers fill the sub vectors from the mama vector
483  for(vecItr=many.begin(),impItr=subImport.begin();
484  vecItr!=many.end();++vecItr,++impItr) {
485  // for ease of access to the destination
486  RCP<Tpetra::MultiVector<ST,LO,GO,NT> > destVec = *vecItr;
487 
488  // extract the map with global indicies from the current vector
489  const Tpetra::Map<LO,GO,NT> & globalMap = *(Teuchos::get_extra_data<RCP<Tpetra::Map<LO,GO,NT> > >(destVec,"globalMap"));
490 
491  // build the import vector as a view on the destination
492  GO destLDA = destVec->getStride();
493  GO destSize = destVec->getGlobalLength()*destVec->getNumVectors();
494  std::vector<ST> destArray(destSize);
495  Teuchos::ArrayView<ST> destVals(destArray);
496  destVec->get1dCopy(destVals,destLDA);
497  Tpetra::MultiVector<ST,LO,GO,NT> importVector(rcpFromRef(globalMap),destVals,destLDA,destVec->getNumVectors());
498 
499  // perform the import
500  importVector.doImport(single,**impItr,Tpetra::INSERT);
501 
502  Tpetra::Import<LO,GO,NT> importer(destVec->getMap(),destVec->getMap());
503  importVector.replaceMap(destVec->getMap());
504  destVec->doImport(importVector,importer,Tpetra::INSERT);
505 
506  }
507 }
508 
509 }
510 } // end namespace Tpetra
511 } // end namespace Teko