Tpetra Matrix/Vector Services Version of the Day
Tpetra_MultiVectorFiller.hpp
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //          Tpetra: Templated Linear Algebra Services Package
00005 //                 Copyright (2008) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 // @HEADER
00041 #ifndef __Tpetra_MultiVectorFiller_hpp
00042 #define __Tpetra_MultiVectorFiller_hpp
00043 
00044 #include <Tpetra_MultiVector.hpp>
00045 #include <Teuchos_CommHelpers.hpp>
00046 #include <iterator>
00047 
00048 namespace {
00049 
00050   // \fn sortAndMergeIn
00051   // \brief Sort and merge newEntries into allEntries, and make unique.
00052   //
00053   // \param allEntries [in/out]: Array in which current entries have
00054   //   already been stored, and in which the new entries are to be
00055   //   stored.  Passed by reference in case we need to resize it.
00056   //
00057   // \param currentEntries [in/out]: Current entries, which we assume
00058   //   have already been sorted and made unique.  Aliases the
00059   //   beginning of allEntries.
00060   //
00061   // \param newEntries [in/out] New entries, which have not yet been
00062   //   sorted or made unique.  This does <i>not</i> alias allEntries.
00063   //
00064   // Sort and make entries of newEntries unique.  Resize allEntries if
00065   // necessary to fit the unique entries of newEntries.  Merge
00066   // newEntries into allEntries and make the results unique.  (This is
00067   // cheaper than sorting the whole array.)
00068   //
00069   // \return A view of all the entries (current and new) in allEntries.
00070   template<class T>
00071   Teuchos::ArrayView<T>
00072   sortAndMergeIn (Teuchos::Array<T>& allEntries,
00073                   Teuchos::ArrayView<T> currentEntries,
00074                   Teuchos::ArrayView<T> newEntries)
00075   {
00076     using Teuchos::ArrayView;
00077     using Teuchos::as;
00078     typedef typename ArrayView<T>::iterator iter_type;
00079     typedef typename ArrayView<T>::size_type size_type;
00080 
00081     std::unique (newEntries.begin(), newEntries.end());
00082     iter_type it = std::unique (newEntries.begin(), newEntries.end());
00083     const size_type numNew = as<size_type> (it - newEntries.begin());
00084     // View of the sorted, made-unique new entries to merge in.
00085     ArrayView<T> newEntriesView = newEntries.view (0, numNew);
00086 
00087     const size_type numCur = currentEntries.size();
00088     if (allEntries.size() < numCur + numNew) {
00089       allEntries.resize (numCur + numNew);
00090     }
00091     ArrayView<T> allView = allEntries.view (0, numCur + numNew);
00092     ArrayView<T> newView = allEntries.view (numCur, numNew); // target of copy
00093 
00094     std::copy (newEntries.begin(), newEntries.end(), newView.begin());
00095     std::inplace_merge (allView.begin(), newView.begin(), allView.end());
00096     iter_type it2 = std::unique (allView.begin(), allView.end());
00097     const size_type numTotal = as<size_type> (it2 - allView.begin());
00098 
00099     return allEntries.view (0, numTotal);
00100   }
00101 
00107   template<class MV>
00108   class MultiVectorFillerData {
00109   public:
00110     typedef typename MV::scalar_type scalar_type;
00111     typedef typename MV::local_ordinal_type local_ordinal_type;
00112     typedef typename MV::global_ordinal_type global_ordinal_type;
00113     typedef typename MV::node_type node_type;
00114 
00115     typedef Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type> map_type;
00116 
00128     MultiVectorFillerData (const Teuchos::RCP<const map_type>& map) :
00129       map_ (map),
00130       numCols_ (0)
00131     {}
00132 
00149     MultiVectorFillerData (const Teuchos::RCP<const map_type>& map,
00150                            const size_t numColumns) :
00151       map_ (map),
00152       numCols_ (numColumns),
00153       sourceIndices_ (numColumns),
00154       sourceValues_ (numColumns)
00155     {}
00156 
00158     void
00159     setNumColumns (const size_t newNumColumns)
00160     {
00161       const size_t oldNumColumns = getNumColumns();
00162       if (newNumColumns >= oldNumColumns) {
00163         for (size_t j = oldNumColumns; j < newNumColumns; ++j) {
00164           sourceIndices_.push_back (Teuchos::Array<global_ordinal_type> ());
00165           sourceValues_.push_back (Teuchos::Array<scalar_type> ());
00166         }
00167       }
00168       else {
00169         // This may not necessarily deallocate any data, but that's OK.
00170         sourceIndices_.resize (newNumColumns);
00171         sourceValues_.resize (newNumColumns);
00172       }
00173       numCols_ = oldNumColumns;
00174     }
00175 
00176     void
00177     sumIntoGlobalValues (Teuchos::ArrayView<const global_ordinal_type> rows,
00178                          size_t column,
00179                          Teuchos::ArrayView<const scalar_type> values)
00180     {
00181       if (column >= getNumColumns()) {
00182         for (size_t j = column; j < getNumColumns(); ++j) {
00183           sourceIndices_.push_back (Teuchos::Array<global_ordinal_type> ());
00184           sourceValues_.push_back (Teuchos::Array<scalar_type> ());
00185         }
00186       }
00187       std::copy (rows.begin(), rows.end(), std::back_inserter (sourceIndices_[column]));
00188       std::copy (values.begin(), values.end(), std::back_inserter (sourceValues_[column]));
00189     }
00190 
00198     void
00199     sumIntoGlobalValues (Teuchos::ArrayView<const global_ordinal_type> rows,
00200                          Teuchos::ArrayView<const scalar_type> values)
00201     {
00202       typedef global_ordinal_type GO;
00203       typedef scalar_type ST;
00204       typedef typename Teuchos::ArrayView<const GO>::const_iterator GoIter;
00205       typedef typename Teuchos::ArrayView<const ST>::const_iterator StIter;
00206 
00207       const size_t numColumns = getNumColumns();
00208       GoIter rowIter = rows.begin();
00209       StIter valIter = values.begin();
00210       for (size_t j = 0; j < numColumns; ++j) {
00211         GoIter rowIterNext = rowIter + numColumns;
00212         StIter valIterNext = valIter + numColumns;
00213         std::copy (rowIter, rowIterNext, std::back_inserter (sourceIndices_[j]));
00214         std::copy (valIter, valIterNext, std::back_inserter (sourceValues_[j]));
00215         rowIter = rowIterNext;
00216         valIter = valIterNext;
00217       }
00218     }
00219 
00244     template<class BinaryFunction>
00245     void
00246     locallyAssemble (MV& X, BinaryFunction& f)
00247     {
00248       using Teuchos::Array;
00249       using Teuchos::ArrayRCP;
00250       using Teuchos::ArrayView;
00251       using Teuchos::RCP;
00252       typedef local_ordinal_type LO;
00253       typedef global_ordinal_type GO;
00254       typedef scalar_type ST;
00255       typedef node_type NT;
00256       typedef Tpetra::Map<LO, GO, NT> map_type;
00257 
00258       RCP<const map_type> map = X.getMap();
00259       Array<LO> localIndices;
00260       const size_t numColumns = getNumColumns();
00261       for (size_t j = 0; j < numColumns; ++j) {
00262         const typename Array<const GO>::size_type numIndices = sourceIndices_[j].size();
00263         // Precompute all the local indices before saving to the
00264         // vector.  Hopefully this gives us a little bit of locality
00265         // in the global->local conversion, at the expense of a little
00266         // more storage.
00267         if (localIndices.size() < numIndices) {
00268           localIndices.resize (numIndices);
00269         }
00270         ArrayView<LO> localIndicesView = localIndices.view (0, numIndices);
00271         ArrayView<const GO> globalIndicesView = sourceIndices_[j].view (0, numIndices);
00272         for (typename ArrayView<const GO>::size_type i = 0; i < numIndices; ++i) {
00273           localIndices[i] = map->getLocalElement (globalIndicesView[i]);
00274         }
00275 
00276         ArrayRCP<ST> X_j = X.getDataNonConst (j);
00277         ArrayView<const ST> localValues = sourceValues_[j].view (0, numIndices);
00278         for (typename ArrayView<const GO>::size_type i = 0; i < numIndices; ++i) {
00279           const LO localInd = localIndices[i];
00280           X_j[localInd] = f (X_j[localInd], localValues[i]);
00281         }
00282       }
00283     }
00284 
00286     void
00287     locallyAssemble (MV& X)
00288     {
00289       std::plus<double> f;
00290       locallyAssemble<std::plus<scalar_type> > (X, f);
00291     }
00292 
00294     void clear() {
00295       Teuchos::Array<Teuchos::Array<global_ordinal_type> > newSourceIndices;
00296       Teuchos::Array<Teuchos::Array<scalar_type> > newSourceValues;
00297       // The standard STL idiom for clearing the contents of a vector completely.
00298       std::swap (sourceIndices_, newSourceIndices);
00299       std::swap (sourceValues_, newSourceValues);
00300     }
00301 
00303     Teuchos::Array<global_ordinal_type>
00304     getSourceIndices () const
00305     {
00306       using Teuchos::Array;
00307       using Teuchos::ArrayView;
00308       using Teuchos::as;
00309       typedef global_ordinal_type GO;
00310       typedef typename Array<GO>::iterator iter_type;
00311       typedef typename Array<GO>::size_type size_type;
00312 
00313       Array<GO> allInds (0); // will resize below
00314       const size_t numCols = getNumColumns();
00315 
00316       if (numCols == 1) {
00317         // Special case for 1 column avoids copying indices twice.
00318         // Pick the size of the array exactly the first time so there
00319         // are at most two allocations (the constructor may choose to
00320         // allocate).
00321         const size_type numNew = sourceIndices_[0].size();
00322         allInds.resize (allInds.size() + numNew);
00323         std::copy (sourceIndices_[0].begin(), sourceIndices_[0].end(),
00324                    allInds.begin());
00325         std::sort (allInds.begin(), allInds.end());
00326         typename Array<GO>::iterator it =
00327           std::unique (allInds.begin(), allInds.end());
00328         const size_type numFinal = as<size_type> (it - allInds.begin());
00329         allInds.resize (numFinal);
00330       }
00331       else {
00332         // Carefully collect all the row indices one column at a time.
00333         // This ensures that the total allocation size in this routine
00334         // is independent of the number of columns.  Also, only sort
00335         // the current column's indices.  Use merges to ensure sorted
00336         // order in the collected final result.
00337         ArrayView<GO> curIndsView = allInds.view (0, 0); // will grow
00338         Array<GO> newInds;
00339         for (size_t j = 0; j < numCols; ++j) {
00340           const size_type numNew = sourceIndices_[j].size();
00341           if (numNew > newInds.size()) {
00342             newInds.resize (numNew);
00343           }
00344           ArrayView<GO> newIndsView = newInds.view (0, numNew);
00345           std::copy (sourceIndices_[j].begin(), sourceIndices_[j].end(),
00346                      newIndsView.begin());
00347           curIndsView = sortAndMergeIn<GO> (allInds, curIndsView, newIndsView);
00348         }
00349       }
00350       return allInds;
00351     }
00352 
00353   private:
00354     Teuchos::RCP<const map_type> map_;
00355     size_t numCols_;
00356     Teuchos::Array<Teuchos::Array<global_ordinal_type> > sourceIndices_;
00357     Teuchos::Array<Teuchos::Array<scalar_type> > sourceValues_;
00358 
00359     size_t getNumColumns() const { return numCols_; }
00360   };
00361 
00367   template<class MV>
00368   class MultiVectorFillerData2 : public Teuchos::Describable {
00369   public:
00370     typedef typename MV::scalar_type scalar_type;
00371     typedef typename MV::local_ordinal_type local_ordinal_type;
00372     typedef typename MV::global_ordinal_type global_ordinal_type;
00373     typedef typename MV::node_type node_type;
00374 
00375     typedef Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type> map_type;
00376 
00387     MultiVectorFillerData2 (const Teuchos::RCP<const map_type>& map,
00388                             const Teuchos::EVerbosityLevel verbLevel=Teuchos::VERB_DEFAULT,
00389                             const Teuchos::RCP<Teuchos::FancyOStream>& out=Teuchos::null) :
00390       map_ (map),
00391       numCols_ (0),
00392       verbLevel_ (verbLevel),
00393       out_ (out)
00394     {}
00395 
00413     MultiVectorFillerData2 (const Teuchos::RCP<const map_type>& map,
00414                             const size_t numColumns,
00415                             const Teuchos::EVerbosityLevel verbLevel=Teuchos::VERB_DEFAULT,
00416                             const Teuchos::RCP<Teuchos::FancyOStream>& out=Teuchos::null) :
00417       map_ (map),
00418       numCols_ (numColumns),
00419       localVec_ (new MV (map, numColumns)),
00420       nonlocalIndices_ (numColumns),
00421       nonlocalValues_ (numColumns),
00422       verbLevel_ (verbLevel),
00423       out_ (out)
00424     {}
00425 
00426     std::string
00427     description() const
00428     {
00429       std::ostringstream oss;
00430       oss << "Tpetra::MultiVectorFillerData2<"
00431           << Teuchos::TypeNameTraits<MV>::name () << ">";
00432       return oss.str();
00433     }
00434 
00435     void
00436     describe (Teuchos::FancyOStream& out,
00437               const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
00438     {
00439       using std::endl;
00440       using Teuchos::Array;
00441       using Teuchos::ArrayRCP;
00442       using Teuchos::ArrayView;
00443       using Teuchos::RCP;
00444       using Teuchos::VERB_DEFAULT;
00445       using Teuchos::VERB_NONE;
00446       using Teuchos::VERB_LOW;
00447       using Teuchos::VERB_MEDIUM;
00448       using Teuchos::VERB_HIGH;
00449       using Teuchos::VERB_EXTREME;
00450 
00451       // Set default verbosity if applicable.
00452       const Teuchos::EVerbosityLevel vl =
00453         (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
00454 
00455       RCP<const Teuchos::Comm<int> > comm = map_->getComm();
00456       const int myImageID = comm->getRank();
00457       const int numImages = comm->getSize();
00458 
00459       if (vl != VERB_NONE) {
00460         // Don't set the tab level unless we're printing something.
00461         Teuchos::OSTab tab (out);
00462 
00463         if (myImageID == 0) { // >= VERB_LOW prints description()
00464           out << this->description() << endl;
00465         }
00466         for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
00467           if (myImageID == imageCtr) {
00468             if (vl != VERB_LOW) {
00469               // At verbosity > VERB_LOW, each process prints something.
00470               out << "Process " << myImageID << ":" << endl;
00471 
00472               Teuchos::OSTab procTab (out);
00473               // >= VERB_MEDIUM: print the local vector length.
00474               out << "local length=" << localVec_->getLocalLength();
00475               if (vl != VERB_MEDIUM) {
00476                 // >= VERB_HIGH: print isConstantStride() and getStride()
00477                 if (localVec_->isConstantStride()) {
00478                   out << ", constant stride=" << localVec_->getStride() << endl;
00479                 }
00480                 else {
00481                   out << ", not constant stride" << endl;
00482                 }
00483                 if (vl == VERB_EXTREME) {
00484                   // VERB_EXTREME: print all the values in the multivector.
00485                   out << "Local values:" << endl;
00486                   ArrayRCP<ArrayRCP<const scalar_type> > X = localVec_->get2dView();
00487                   for (size_t i = 0; i < localVec_->getLocalLength(); ++i) {
00488                     for (size_t j = 0; j < localVec_->getNumVectors(); ++j) {
00489                       out << X[j][i];
00490                       if (j < localVec_->getNumVectors() - 1) {
00491                         out << " ";
00492                       }
00493                     } // for each column
00494                     out << endl;
00495                   } // for each row
00496 
00497                   out << "Nonlocal indices and values:" << endl;
00498                   for (size_t j = 0; j < (size_t)nonlocalIndices_.size(); ++j) {
00499                     ArrayView<const global_ordinal_type> inds = nonlocalIndices_[j]();
00500                     ArrayView<const scalar_type> vals = nonlocalValues_[j]();
00501 
00502                     for (typename ArrayView<const global_ordinal_type>::size_type k = 0; k < inds.size(); ++k) {
00503                       out << "X(" << inds[k] << "," << j << ") = " << vals[k] << endl;
00504                     }
00505                   }
00506                 } // if vl == VERB_EXTREME
00507               } // if (vl != VERB_MEDIUM)
00508               else { // vl == VERB_LOW
00509                 out << endl;
00510               }
00511             } // if vl != VERB_LOW
00512           } // if it is my process' turn to print
00513         } // for each process in the communicator
00514       } // if vl != VERB_NONE
00515     }
00516 
00518     Teuchos::Array<global_ordinal_type>
00519     getSourceIndices () const
00520     {
00521       using Teuchos::Array;
00522       using Teuchos::ArrayView;
00523       using Teuchos::RCP;
00524       using Teuchos::rcp;
00525       using Tpetra::global_size_t;
00526       typedef global_ordinal_type GO;
00527 
00528       // Get the nonlocal row indices, sorted and made unique.
00529       // It's fair to assume that these are not contiguous.
00530       Array<GO> nonlocalIndices = getSortedUniqueNonlocalIndices();
00531 
00532       // Get the local row indices, not necessarily sorted or unique.
00533       ArrayView<const GO> localIndices = getLocalIndices ();
00534 
00535       // Copy the local indices into the full indices array, and sort
00536       // them there.  We'll merge in the nonlocal indices below.  This
00537       // can be more efficient than just sorting all the indices, if
00538       // there are a lot of nonlocal indices.
00539       Array<GO> indices (localIndices.size() + nonlocalIndices.size());
00540       ArrayView<GO> localIndView = indices.view (0, localIndices.size());
00541       std::copy (localIndices.begin(), localIndices.end(), localIndView.begin());
00542       std::sort (localIndView.begin(), localIndView.end());
00543 
00544       // Merge the local and nonlocal indices.
00545       if (nonlocalIndices.size() > 0) {
00546         typedef typename ArrayView<GO>::iterator iter_type;
00547 
00548         iter_type middle = localIndView.end();
00549         // We need a view, because std::inplace_merge needs all its
00550         // iterator inputs to have the same type.  Debug mode builds
00551         // are pickier than release mode builds, because the iterators
00552         // in a debug mode build are of a different type that does
00553         // run-time checking (they aren't just raw pointers).
00554         ArrayView<GO> indView = indices.view (0, indices.size());
00555 
00556         //iter_type middle = &indices[nonlocalIndices.size()];
00557         std::copy (nonlocalIndices.begin(), nonlocalIndices.end(), middle);
00558         std::inplace_merge (indView.begin(), middle, indView.end());
00559       }
00560       return indices;
00561     }
00562 
00572     void
00573     setNumColumns (const size_t newNumColumns)
00574     {
00575       using Teuchos::Array;
00576       using Teuchos::Range1D;
00577       using Teuchos::RCP;
00578       typedef global_ordinal_type GO;
00579       typedef scalar_type ST;
00580 
00581       const size_t oldNumColumns = numCols_;
00582       if (newNumColumns == oldNumColumns) {
00583         return; // No side effects if no change.
00584       }
00585 
00586       RCP<MV> newLocalVec;
00587       if (newNumColumns > oldNumColumns) {
00588         newLocalVec = rcp (new MV (map_, newNumColumns));
00589         // Assign the contents of the old local multivector to the
00590         // first oldNumColumns columns of the new local multivector,
00591         // then get rid of the old local multivector.
00592         RCP<MV> newLocalVecView =
00593           newLocalVec->subViewNonConst (Range1D (0, oldNumColumns-1));
00594         *newLocalVecView = *localVec_;
00595       }
00596       else {
00597         if (newNumColumns == 0) {
00598           // Tpetra::MultiVector doesn't let you construct a
00599           // multivector with zero columns.
00600           newLocalVec = Teuchos::null;
00601         }
00602         else {
00603           newLocalVec =
00604             localVec_->subViewNonConst (Range1D (0, newNumColumns-1));
00605         }
00606       }
00607 
00608       // Leave most side effects until the end, for exception safety.
00609       nonlocalIndices_.resize (newNumColumns);
00610       nonlocalValues_.resize (newNumColumns);
00611       localVec_ = newLocalVec;
00612       numCols_ = newNumColumns;
00613     }
00614 
00616     void
00617     sumIntoGlobalValues (Teuchos::ArrayView<const global_ordinal_type> rows,
00618                          size_t columnIndex,
00619                          Teuchos::ArrayView<const scalar_type> values)
00620     {
00621       using Teuchos::ArrayRCP;
00622       using Teuchos::ArrayView;
00623       typedef local_ordinal_type LO;
00624       typedef global_ordinal_type GO;
00625       typedef scalar_type ST;
00626 
00627       if (columnIndex >= getNumColumns()) {
00628         // Automatically expand the number of columns.  This
00629         // implicitly ensures that localVec_ is not null.
00630         setNumColumns (columnIndex + 1);
00631       }
00632 
00633       typename ArrayView<const GO>::const_iterator rowIter = rows.begin();
00634       typename ArrayView<const ST>::const_iterator valIter = values.begin();
00635       for ( ; rowIter != rows.end() && valIter != values.end(); ++rowIter, ++valIter) {
00636         const GO globalRowIndex = *rowIter;
00637         // Converting from global to local index could be logarithmic
00638         // in the number of global indices that this process owns,
00639         // depending on the Map implementation.  However, the lookup
00640         // allows us to store data in the local multivector, rather
00641         // than in a separate data structure.
00642         const LO localRowIndex = map_->getLocalElement (globalRowIndex);
00643         if (localRowIndex == Teuchos::OrdinalTraits<LO>::invalid()) {
00644           nonlocalIndices_[columnIndex].push_back (globalRowIndex);
00645           nonlocalValues_[columnIndex].push_back (*valIter);
00646         }
00647         else {
00648           // FIXME (mfh 27 Feb 2012) This will be very slow for GPU
00649           // Node types.  In that case, we should hold on to the view
00650           // of localVec_ as long as the number of columns doesn't
00651           // change, and make modifications to the view until
00652           // localAssemble() is called.
00653           ArrayRCP<ST> X_j = localVec_->getDataNonConst (columnIndex);
00654           // FIXME (mfh 27 Feb 2012) Allow different combine modes.
00655           // The current combine mode just adds to the current value
00656           // at that location.
00657           X_j[localRowIndex] += *valIter;
00658         }
00659       }
00660     }
00661 
00671     void
00672     sumIntoGlobalValues (Teuchos::ArrayView<const global_ordinal_type> rows,
00673                          Teuchos::ArrayView<const scalar_type> values)
00674     {
00675       using Teuchos::ArrayView;
00676       typedef typename ArrayView<const global_ordinal_type>::size_type size_type;
00677 
00678       const size_t numCols = getNumColumns();
00679       for (size_t j = 0; j < numCols; ++j) {
00680         const size_type offset = numCols*j;
00681         const size_type len = numCols;
00682         sumIntoGlobalValues (rows.view (offset, len), j, values.view (offset, len));
00683       }
00684     }
00685 
00710     template<class BinaryFunction>
00711     void
00712     locallyAssemble (MV& X, BinaryFunction& f)
00713     {
00714       using Teuchos::ArrayRCP;
00715       using Teuchos::ArrayView;
00716       using Teuchos::FancyOStream;
00717       using Teuchos::getFancyOStream;
00718       using Teuchos::oblackholestream;
00719       using Teuchos::RCP;
00720       using Teuchos::rcp;
00721       using std::endl;
00722 
00723       typedef local_ordinal_type LO;
00724       typedef global_ordinal_type GO;
00725       typedef scalar_type ST;
00726 
00727       // Default output stream prints nothing.
00728       RCP<FancyOStream> out = out_.is_null() ?
00729         getFancyOStream (rcp (new oblackholestream)) : out_;
00730 
00731       Teuchos::OSTab tab (out);
00732       *out << "locallyAssemble:" << endl;
00733 
00734       RCP<const map_type> srcMap = X.getMap();
00735       ArrayView<const GO> localIndices = map_->getNodeElementList ();
00736 
00737       for (size_t j = 0; j < X.getNumVectors(); ++j) {
00738         ArrayRCP<ST> X_j = X.getDataNonConst (j);
00739 
00740         // First add all the local data into X_j.
00741         ArrayRCP<const ST> local_j = localVec_->getDataNonConst (j);
00742         for (typename ArrayView<const GO>::const_iterator it = localIndices.begin();
00743              it != localIndices.end(); ++it) {
00744           const LO rowIndLocal = map_->getLocalElement (*it);
00745           const LO rowIndX = srcMap->getLocalElement (*it);
00746 
00747           TEUCHOS_TEST_FOR_EXCEPTION(rowIndX == Teuchos::OrdinalTraits<LO>::invalid(),
00748             std::invalid_argument, "locallyAssemble(): Input multivector X does "
00749             "not own the global index " << *it << ".  This probably means that "
00750             "X was not constructed with the right Map.");
00751           X_j[rowIndX] = f (X_j[rowIndX], local_j[rowIndLocal]);
00752         }
00753 
00754         // Now add the nonlocal data into X_j.
00755         ArrayView<const GO> nonlocalIndices = nonlocalIndices_[j]();
00756         typename ArrayView<const GO>::const_iterator indexIter = nonlocalIndices.begin();
00757         ArrayView<const ST> nonlocalValues = nonlocalValues_[j]();
00758         typename ArrayView<const ST>::const_iterator valueIter = nonlocalValues.begin();
00759         for ( ; indexIter != nonlocalIndices.end() && valueIter != nonlocalValues.end();
00760               ++indexIter, ++valueIter) {
00761           const LO rowIndX = srcMap->getLocalElement (*indexIter);
00762           X_j[rowIndX] = f (X_j[rowIndX], *valueIter);
00763         }
00764       }
00765 
00766       *out << "Locally assembled vector:" << endl;
00767       X.describe (*out, Teuchos::VERB_EXTREME);
00768     }
00769 
00771     void
00772     locallyAssemble (MV& X)
00773     {
00774       std::plus<double> f;
00775       locallyAssemble<std::plus<scalar_type> > (X, f);
00776     }
00777 
00783     void clear() {
00784       Teuchos::Array<Teuchos::Array<global_ordinal_type> > newNonlocalIndices;
00785       Teuchos::Array<Teuchos::Array<scalar_type> > newNonlocalValues;
00786       // The standard STL idiom for clearing the contents of a vector
00787       // completely.  Setting the size to zero may not necessarily
00788       // deallocate data.
00789       std::swap (nonlocalIndices_, newNonlocalIndices);
00790       std::swap (nonlocalValues_, newNonlocalValues);
00791 
00792       // Don't actually deallocate the multivector of local entries.
00793       // Just fill it with zero.  This is because the caller hasn't
00794       // reset the number of columns.
00795       if (! localVec_.is_null()) {
00796         localVec_->putScalar (Teuchos::ScalarTraits<scalar_type>::zero());
00797       }
00798     }
00799 
00800   private:
00802     Teuchos::RCP<const map_type> map_;
00803 
00805     size_t numCols_;
00806 
00808     Teuchos::RCP<MV> localVec_;
00809 
00811     Teuchos::Array<Teuchos::Array<global_ordinal_type> > nonlocalIndices_;
00812 
00814     Teuchos::Array<Teuchos::Array<scalar_type> > nonlocalValues_;
00815 
00817     Teuchos::EVerbosityLevel verbLevel_;
00818 
00820     Teuchos::RCP<Teuchos::FancyOStream> out_;
00821 
00823     size_t getNumColumns() const { return numCols_; }
00824 
00826     Teuchos::ArrayView<const global_ordinal_type>
00827     getLocalIndices() const
00828     {
00829       return map_->getNodeElementList ();
00830     }
00831 
00833     Teuchos::Array<global_ordinal_type>
00834     getSortedUniqueNonlocalIndices() const
00835     {
00836       using Teuchos::Array;
00837       using Teuchos::ArrayView;
00838       using Teuchos::as;
00839       typedef global_ordinal_type GO;
00840       typedef typename Array<GO>::size_type size_type;
00841 
00842       Array<GO> allNonlocals (0); // will resize below
00843       const size_t numCols = getNumColumns();
00844 
00845       if (numCols == 1) {
00846         // Special case for 1 column avoids copying indices twice.
00847         // Pick the size of the array exactly the first time so there
00848         // are at most two allocations (the constructor may choose to
00849         // allocate).
00850         const size_type numNew = nonlocalIndices_[0].size();
00851         allNonlocals.resize (allNonlocals.size() + numNew);
00852         std::copy (nonlocalIndices_[0].begin(), nonlocalIndices_[0].end(),
00853                    allNonlocals.begin());
00854         std::sort (allNonlocals.begin(), allNonlocals.end());
00855         typename Array<GO>::iterator it =
00856           std::unique (allNonlocals.begin(), allNonlocals.end());
00857         const size_type numFinal = as<size_type> (it - allNonlocals.begin());
00858         allNonlocals.resize (numFinal);
00859       }
00860       else {
00861         // Carefully collect all the row indices one column at a time.
00862         // This ensures that the total allocation size in this routine
00863         // is independent of the number of columns.  Also, only sort
00864         // the current column's indices.  Use merges to ensure sorted
00865         // order in the collected final result.
00866         ArrayView<GO> curNonlocalsView = allNonlocals.view (0, 0); // will grow
00867         Array<GO> newNonlocals;
00868         for (size_t j = 0; j < numCols; ++j) {
00869           const size_type numNew = nonlocalIndices_[j].size();
00870           if (numNew > newNonlocals.size()) {
00871             newNonlocals.resize (numNew);
00872           }
00873           ArrayView<GO> newNonlocalsView = newNonlocals.view (0, numNew);
00874           std::copy (nonlocalIndices_[j].begin(), nonlocalIndices_[j].end(),
00875                      newNonlocalsView.begin());
00876           curNonlocalsView = sortAndMergeIn<GO> (allNonlocals, curNonlocalsView,
00877                                                  newNonlocalsView);
00878         }
00879       }
00880       return allNonlocals;
00881     }
00882   };
00883 } // namespace (anonymous)
00884 
00885 namespace Tpetra {
00886 
00911   template<class MV>
00912   class MultiVectorFiller {
00913   public:
00914     typedef typename MV::scalar_type scalar_type;
00915     typedef typename MV::local_ordinal_type local_ordinal_type;
00916     typedef typename MV::global_ordinal_type global_ordinal_type;
00917     typedef typename MV::node_type node_type;
00918     typedef Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type> map_type;
00919 
00949     MultiVectorFiller (const Teuchos::RCP<const map_type>& map,
00950                        const size_t numCols);
00951 
00971     void globalAssemble (MV& X_out, const bool forceReuseMap = false);
00972 
00973     void
00974     describe (Teuchos::FancyOStream& out,
00975               const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
00976     {
00977       data_.describe (out, verbLevel);
00978     }
00979 
00990     void
00991     sumIntoGlobalValues (Teuchos::ArrayView<const global_ordinal_type> rows,
00992                          size_t column,
00993                          Teuchos::ArrayView<const scalar_type> values)
00994     {
00995       data_.sumIntoGlobalValues (rows, column, values);
00996     }
00997 
01013     void
01014     sumIntoGlobalValues (Teuchos::ArrayView<const global_ordinal_type> rows,
01015                          Teuchos::ArrayView<const scalar_type> values)
01016     {
01017       data_.sumIntoGlobalValues (rows, values);
01018     }
01019 
01020   private:
01022     Teuchos::RCP<const map_type> ctorMap_;
01023 
01028     Teuchos::RCP<const map_type> sourceMap_;
01029 
01035     Teuchos::RCP<const map_type> targetMap_;
01036 
01049     Teuchos::RCP<MV> sourceVec_;
01050 
01055     MultiVectorFillerData2<MV> data_;
01056 
01057     typedef Tpetra::Export<local_ordinal_type, global_ordinal_type, node_type> export_type;
01058     Teuchos::RCP<export_type> exporter_;
01059 
01065     void locallyAssemble (MV& X_in) {
01066       data_.locallyAssemble (X_in);
01067     }
01068 
01070     Teuchos::Array<global_ordinal_type> getSourceIndices () const {
01071       return data_.getSourceIndices();
01072     }
01073 
01082     Teuchos::RCP<const map_type>
01083     computeSourceMap (const global_ordinal_type indexBase,
01084                       const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
01085                       const Teuchos::RCP<node_type>& node);
01086   };
01087 
01088   template<class MV>
01089   MultiVectorFiller<MV>::MultiVectorFiller (const Teuchos::RCP<const typename MultiVectorFiller<MV>::map_type>& map, const size_t numCols)
01090     : ctorMap_ (map),
01091       sourceMap_ (Teuchos::null),
01092       targetMap_ (Teuchos::null),
01093       data_ (map, numCols),
01094       exporter_ (Teuchos::null)
01095   {}
01096 
01097   template<class MV>
01098   Teuchos::RCP<const typename MultiVectorFiller<MV>::map_type>
01099   MultiVectorFiller<MV>::
01100   computeSourceMap (const global_ordinal_type indexBase,
01101                     const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
01102                     const Teuchos::RCP<node_type>& node)
01103   {
01104     using Teuchos::Array;
01105     using Teuchos::ArrayView;
01106     using Teuchos::rcp;
01107     typedef global_ordinal_type GO;
01108 
01109     Array<GO> indices = getSourceIndices ();
01110 
01111     // Passing "invalid" for the numGlobalElements argument of the Map
01112     // constructor tells the Map to compute the global number of
01113     // elements itself.
01114     const Tpetra::global_size_t invalid =
01115       Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
01116     return rcp (new map_type (invalid, indices, indexBase, comm, node));
01117   }
01118 
01119   template<class MV>
01120   void
01121   MultiVectorFiller<MV>::globalAssemble (MV& X_out, const bool forceReuseMap)
01122   {
01123     using Teuchos::ArrayView;
01124     using Teuchos::Array;
01125     using Teuchos::Range1D;
01126     using Teuchos::RCP;
01127     using Teuchos::rcp;
01128     typedef global_ordinal_type GO;
01129 
01130     const size_t numVecs = X_out.getNumVectors();
01131 
01132     if (numVecs == 0) {
01133       // Nothing to do!  Of course, this does not check for whether
01134       // X_out has the right number of rows.  That's OK, though.
01135       // Compare to the definition of the BLAS' _AXPY for an input
01136       // vector containing NaNs, but multiplied by alpha=0.
01137       return;
01138     }
01139     //
01140     // Get the target Map of the Export.  If X_out's Map is the same
01141     // as the target Map from last time, then we may be able to
01142     // recycle the Export from last time, if the source Map is also
01143     // the same.
01144     //
01145     RCP<const map_type> targetMap;
01146     bool assumeSameTargetMap = false;
01147     if (targetMap_.is_null()) {
01148       targetMap_ = X_out.getMap();
01149       targetMap = targetMap_;
01150       assumeSameTargetMap = false;
01151     }
01152     else {
01153       if (forceReuseMap) {
01154         targetMap = targetMap_;
01155         assumeSameTargetMap = true;
01156       }
01157       else  {
01158         // If X_out's Map is the same as targetMap_, we may be able to
01159         // reuse the Export.  Constructing the Export may be more
01160         // expensive than calling isSameAs() (which requires just a
01161         // few reductions and reading through the lists of owned
01162         // global indices), so it's worth checking.
01163         if (targetMap_->isSameAs (*(X_out.getMap()))) {
01164           assumeSameTargetMap = true;
01165           targetMap = targetMap_;
01166         }
01167       }
01168     }
01169     //
01170     // Get the source Map of the Export.  If the source Map of the
01171     // Export is the same as last time, then we may be able to recycle
01172     // the Export from last time, if the target Map is also the same.
01173     //
01174     RCP<const map_type> sourceMap;
01175     bool computedSourceMap = false;
01176     {
01177       if (forceReuseMap && ! sourceMap_.is_null()) {
01178         sourceMap = sourceMap_;
01179       }
01180       else {
01181         sourceMap = computeSourceMap (ctorMap_->getIndexBase(),
01182                                       ctorMap_->getComm(),
01183                                       ctorMap_->getNode());
01184         computedSourceMap = true;
01185       }
01186     }
01187     if (computedSourceMap) {
01188       sourceMap_ = sourceMap;
01189     }
01190     //
01191     // Now that we have the source and target Maps of the Export, we
01192     // can check whether we can recycle the Export from last time.
01193     //
01194     const bool mustComputeExport =
01195       (exporter_.is_null() || (assumeSameTargetMap && ! computedSourceMap));
01196     if (mustComputeExport) {
01197       exporter_ = rcp (new export_type (sourceMap_, targetMap_));
01198     }
01199 
01200     // Source multivector for the Export.
01201     RCP<MV> X_in;
01202     const bool mustReallocateVec = sourceVec_.is_null() ||
01203       sourceVec_->getNumVectors() < numVecs || ! assumeSameTargetMap;
01204 
01205     if (mustReallocateVec) {
01206       // Reallocating nonlocalVec_ ensures that it has the right Map.
01207       sourceVec_ = rcp (new MV (sourceMap_, numVecs));
01208       X_in = sourceVec_;
01209     } else {
01210       if (sourceVec_->getNumVectors() == numVecs) {
01211         X_in = sourceVec_;
01212       } else { // sourceVec_ has more vectors than needed.
01213         X_in = sourceVec_->subViewNonConst (Range1D (0, numVecs-1));
01214       }
01215     }
01216 
01217     // "Locally assemble" the data into X_in by summing together
01218     // entries with the same indices.
01219     locallyAssemble (*X_in);
01220 
01221     // Do the Export.
01222     const Tpetra::CombineMode combineMode = Tpetra::ADD;
01223     X_out.doExport (*X_in, *exporter_, combineMode);
01224   }
01225 
01226   namespace Test {
01227 
01233     template<class MV>
01234     class MultiVectorFillerTester {
01235     public:
01236       typedef typename MV::scalar_type scalar_type;
01237       typedef typename MV::local_ordinal_type local_ordinal_type;
01238       typedef typename MV::global_ordinal_type global_ordinal_type;
01239       typedef typename MV::node_type node_type;
01240       typedef Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type> map_type;
01241 
01250       static void
01251       testSameMap (const Teuchos::RCP<const map_type>& targetMap,
01252                    const global_ordinal_type eltSize, // Must be odd
01253                    const size_t numCols,
01254                    const Teuchos::RCP<Teuchos::FancyOStream>& outStream=Teuchos::null,
01255                    const Teuchos::EVerbosityLevel verbosityLevel=Teuchos::VERB_DEFAULT)
01256       {
01257         using Teuchos::Array;
01258         using Teuchos::ArrayRCP;
01259         using Teuchos::ArrayView;
01260         using Teuchos::as;
01261         using Teuchos::Comm;
01262         using Teuchos::FancyOStream;
01263         using Teuchos::getFancyOStream;
01264         using Teuchos::oblackholestream;
01265         using Teuchos::ptr;
01266         using Teuchos::RCP;
01267         using Teuchos::rcp;
01268         using Teuchos::rcpFromRef;
01269         using Teuchos::REDUCE_SUM;
01270         using Teuchos::reduceAll;
01271         using std::cerr;
01272         using std::endl;
01273 
01274         typedef local_ordinal_type LO;
01275         typedef global_ordinal_type GO;
01276         typedef scalar_type ST;
01277         typedef Teuchos::ScalarTraits<ST> STS;
01278 
01279         TEUCHOS_TEST_FOR_EXCEPTION(eltSize % 2 == 0, std::invalid_argument,
01280           "Element size (eltSize) argument must be odd.");
01281         TEUCHOS_TEST_FOR_EXCEPTION(numCols == 0, std::invalid_argument,
01282           "Number of columns (numCols) argument must be nonzero.");
01283         // Default behavior is to print nothing out.
01284         RCP<FancyOStream> out = outStream.is_null() ?
01285           getFancyOStream (rcp (new oblackholestream)) : outStream;
01286         const Teuchos::EVerbosityLevel verbLevel =
01287           (verbosityLevel == Teuchos::VERB_DEFAULT) ?
01288           Teuchos::VERB_NONE : verbosityLevel;
01289 
01290         //RCP<MV> X = rcp (new MV (targetMap, numCols));
01291 
01292         const GO indexBase = targetMap->getIndexBase();
01293         Array<GO> rows (eltSize);
01294         Array<ST> values (eltSize);
01295         std::fill (values.begin(), values.end(), STS::one());
01296 
01297         // Make this a pointer so we can free its contents, in case
01298         // those contents depend on the input to globalAssemble().
01299         RCP<MultiVectorFiller<MV> > filler =
01300           rcp (new MultiVectorFiller<MV> (targetMap, numCols));
01301 
01302         TEUCHOS_TEST_FOR_EXCEPTION(! targetMap->isContiguous(),
01303           std::invalid_argument, "MultiVectorFiller test currently only works "
01304           "for contiguous Maps.");
01305 
01306         const GO minGlobalIndex = targetMap->getMinGlobalIndex();
01307         const GO maxGlobalIndex = targetMap->getMaxGlobalIndex();
01308         const GO minAllGlobalIndex = targetMap->getMinAllGlobalIndex();
01309         const GO maxAllGlobalIndex = targetMap->getMaxAllGlobalIndex();
01310         for (size_t j = 0; j < numCols; ++j) {
01311           for (GO i = minGlobalIndex; i <= maxGlobalIndex; ++i) {
01312             // Overlap over processes, without running out of bounds.
01313             const GO start = std::max (i - eltSize/2, minAllGlobalIndex);
01314             const GO end = std::min (i + eltSize/2, maxAllGlobalIndex);
01315             const GO len = end - start + 1;
01316 
01317             TEUCHOS_TEST_FOR_EXCEPTION(len > eltSize, std::logic_error,
01318               "At start,end = " << start << "," << end << ", len = " << len
01319               << " > eltSize = " << eltSize << ".");
01320 
01321             for (GO k = 0; k < len; ++k) {
01322               rows[k] = start + k;
01323             }
01324             if (verbLevel == Teuchos::VERB_EXTREME) {
01325               *out << "Inserting: "
01326                    << Teuchos::toString (rows.view(0,len)) << ", "
01327                    << Teuchos::toString (values.view(0, len)) << std::endl;
01328             }
01329             filler->sumIntoGlobalValues (rows.view(0, len), j, values.view(0, len));
01330           }
01331         }
01332 
01333         if (verbLevel == Teuchos::VERB_EXTREME) {
01334           *out << "Filler:" << std::endl;
01335           filler->describe (*out, verbLevel);
01336           *out << std::endl;
01337         }
01338 
01339         MV X_out (targetMap, numCols);
01340         filler->globalAssemble (X_out);
01341         filler = Teuchos::null;
01342 
01343         if (verbLevel == Teuchos::VERB_EXTREME) {
01344           *out << "X_out:" << std::endl;
01345           X_out.describe (*out, verbLevel);
01346           *out << std::endl;
01347         }
01348 
01349         // Create multivector for comparison.
01350         MV X_expected (targetMap, numCols);
01351         const scalar_type two = STS::one() + STS::one();
01352         for (size_t j = 0; j < numCols; ++j) {
01353           ArrayRCP<ST> X_j = X_expected.getDataNonConst (j);
01354 
01355           // Each entry of the column should have the value eltSize,
01356           // except for the first and last few entries in the whole
01357           // column (globally, not locally).
01358           for (GO i = minGlobalIndex; i <= maxGlobalIndex; ++i) {
01359             const LO localIndex = targetMap->getLocalElement (i);
01360             TEUCHOS_TEST_FOR_EXCEPTION(i == Teuchos::OrdinalTraits<LO>::invalid(),
01361               std::logic_error, "Global index " << i << " is not in the "
01362               "multivector's Map.");
01363 
01364             if (i <= minAllGlobalIndex + eltSize/2) {
01365               X_j[localIndex] = two + as<ST>(i) - as<ST>(indexBase);
01366             }
01367             else if (i >= maxAllGlobalIndex - eltSize/2) {
01368               X_j[localIndex] = two + as<ST>(maxAllGlobalIndex) - as<ST>(i);
01369             }
01370             else {
01371               X_j[localIndex] = as<ST>(eltSize);
01372             }
01373           } // for each global index which my process owns
01374         } // for each column of the multivector
01375 
01376         if (verbLevel == Teuchos::VERB_EXTREME) {
01377           *out << "X_expected:" << std::endl;
01378           X_expected.describe (*out, verbLevel);
01379           *out << std::endl;
01380         }
01381 
01382         Array<GO> errorLocations;
01383         for (size_t j = 0; j < numCols; ++j) {
01384           ArrayRCP<const ST> X_out_j = X_out.getData (j);
01385           ArrayRCP<const ST> X_expected_j = X_expected.getData (j);
01386 
01387           // Each entry of the column should have the value eltSize,
01388           // except for the first and last few entries in the whole
01389           // column (globally, not locally).
01390           for (GO i = minGlobalIndex; i <= maxGlobalIndex; ++i) {
01391             const LO localIndex = targetMap->getLocalElement (i);
01392             TEUCHOS_TEST_FOR_EXCEPTION(i == Teuchos::OrdinalTraits<LO>::invalid(),
01393               std::logic_error, "Global index " << i << " is not in the "
01394               "multivector's Map.");
01395 
01396             // The floating-point additions should be exact in this
01397             // case, except for very large values of eltSize.
01398             if (X_out_j[localIndex] != X_expected_j[localIndex]) {
01399               errorLocations.push_back (i);
01400             }
01401           } // for each global index which my process owns
01402 
01403           const typename Array<GO>::size_type localNumErrors = errorLocations.size();
01404           typename Array<GO>::size_type globalNumErrors = 0;
01405           RCP<const Comm<int> > comm = targetMap->getComm();
01406           reduceAll (*comm, REDUCE_SUM, localNumErrors, ptr (&globalNumErrors));
01407 
01408           if (globalNumErrors != 0) {
01409             std::ostringstream os;
01410             os << "Proc " << comm->getRank() << ": " << localNumErrors
01411                << " incorrect value" << (localNumErrors != 1 ? "s" : "")
01412                << ".  Error locations: [ ";
01413             std::copy (errorLocations.begin(), errorLocations.end(),
01414                        std::ostream_iterator<GO> (os, " "));
01415             os << "].";
01416             // Iterate through all processes in the communicator,
01417             // printing out each process' local errors.
01418             for (int p = 0; p < comm->getSize(); ++p) {
01419               if (p == comm->getRank()) {
01420                 cerr << os.str() << endl;
01421               }
01422               // Barriers to let output finish.
01423               comm->barrier();
01424               comm->barrier();
01425               comm->barrier();
01426             }
01427             TEUCHOS_TEST_FOR_EXCEPTION(globalNumErrors != 0, std::logic_error,
01428               "Over all procs: " << globalNumErrors << " total error"
01429               << (globalNumErrors != 1 ? "s" : "") << ".");
01430           } // if there were any errors in column j
01431         } // for each column j
01432       }
01433     };
01434 
01436     template<class ScalarType,
01437              class LocalOrdinalType,
01438              class GlobalOrdinalType,
01439              class NodeType>
01440     void
01441     testMultiVectorFiller (const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
01442                            const Teuchos::RCP<NodeType>& node,
01443                            const size_t unknownsPerNode,
01444                            const GlobalOrdinalType unknownsPerElt,
01445                            const size_t numCols,
01446                            const Teuchos::RCP<Teuchos::FancyOStream>& outStream,
01447                            const Teuchos::EVerbosityLevel verbLevel)
01448     {
01449       using Tpetra::createContigMapWithNode;
01450       using Teuchos::FancyOStream;
01451       using Teuchos::getFancyOStream;
01452       using Teuchos::oblackholestream;
01453       using Teuchos::ParameterList;
01454       using Teuchos::RCP;
01455       using Teuchos::rcp;
01456       using std::cerr;
01457       using std::endl;
01458 
01459       typedef ScalarType ST;
01460       typedef LocalOrdinalType LO;
01461       typedef GlobalOrdinalType GO;
01462       typedef NodeType NT;
01463       typedef Tpetra::Map<LO, GO, NT> MT;
01464       typedef Tpetra::MultiVector<ST, LO, GO, NT> MV;
01465 
01466       RCP<FancyOStream> out = outStream.is_null() ?
01467         getFancyOStream (rcp (new oblackholestream)) : outStream;
01468       const Tpetra::global_size_t invalid =
01469         Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
01470       RCP<const MT> targetMap =
01471         createContigMapWithNode<LO, GO, NT> (invalid, unknownsPerNode, comm, node);
01472 
01473       std::ostringstream os;
01474       int success = 1;
01475       try {
01476         MultiVectorFillerTester<MV>::testSameMap (targetMap, unknownsPerElt, numCols, out, verbLevel);
01477         success = 1;
01478       } catch (std::exception& e) {
01479         success = 0;
01480         os << e.what();
01481       }
01482 
01483       for (int p = 0; p < comm->getSize(); ++p) {
01484         if (p == comm->getRank()) {
01485           *out << "On process " << comm->getRank() << ": " << os.str() << endl;
01486         }
01487         comm->barrier();
01488         comm->barrier();
01489         comm->barrier();
01490       }
01491     }
01492 
01498     void
01499     testSortAndMergeIn ();
01500 
01501   } // namespace Test
01502 } // namespace Tpetra
01503 
01504 
01505 #endif // __Tpetra_MultiVectorFiller_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines