Tpetra Matrix/Vector Services Version of the Day
Tpetra_Vector_def.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 
00042 #ifndef TPETRA_VECTOR_DEF_HPP
00043 #define TPETRA_VECTOR_DEF_HPP
00044 
00045 #include <Kokkos_DefaultArithmetic.hpp>
00046 #include <Kokkos_NodeTrace.hpp>
00047 #include <Tpetra_MultiVector.hpp>
00048 
00049 #ifdef DOXYGEN_USE_ONLY
00050   #include "Tpetra_Vector_decl.hpp"
00051 #endif
00052 
00053 namespace Tpetra {
00054 
00055   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00056   Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Vector(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map, bool zeroOut)
00057     : MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map,1,zeroOut) {
00058   }
00059 
00060   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00061   Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Vector(const Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &source)
00062   : MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(source) {
00063   }
00064 
00065   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00066   Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Vector(
00067                               const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map,
00068                               const ArrayRCP<Scalar> &view, EPrivateHostViewConstructor /* dummy */)
00069   : MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map,view,view.size(),1,HOST_VIEW_CONSTRUCTOR) {
00070   }
00071 
00072   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00073   Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Vector(
00074                               const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map,
00075                               const ArrayRCP<Scalar> &view, EPrivateComputeViewConstructor /* dummy */)
00076   : MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map,view,view.size(),1,COMPUTE_VIEW_CONSTRUCTOR) {
00077   }
00078 
00079   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00080   Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Vector(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map, const ArrayView<const Scalar> &values)
00081   : MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map,values,values.size(),1) {
00082   }
00083 
00084   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00085   Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::~Vector() {}
00086 
00087   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00088   void Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::replaceGlobalValue(GlobalOrdinal globalRow, const Scalar &value) {
00089     this->MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::replaceGlobalValue(globalRow,0,value);
00090   }
00091 
00092   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00093   void Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::sumIntoGlobalValue(GlobalOrdinal globalRow, const Scalar &value) {
00094     this->MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::sumIntoGlobalValue(globalRow,0,value);
00095   }
00096 
00097   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00098   void Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::replaceLocalValue(LocalOrdinal myRow, const Scalar &value) {
00099     this->MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::replaceLocalValue(myRow,0,value);
00100   }
00101 
00102   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00103   void Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::sumIntoLocalValue(LocalOrdinal myRow, const Scalar &value) {
00104     this->MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::sumIntoLocalValue(myRow,0,value);
00105   }
00106 
00107   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00108   void Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::get1dCopy(ArrayView<Scalar> A) const {
00109     size_t lda = this->getLocalLength();
00110     this->get1dCopy(A,lda);
00111   }
00112 
00113   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00114   Scalar Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::dot(const Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &a) const {
00115     using Teuchos::outArg;
00116 #ifdef HAVE_TPETRA_DEBUG
00117     TEUCHOS_TEST_FOR_EXCEPTION( !this->getMap()->isCompatible(*a.getMap()), std::runtime_error,
00118         "Tpetra::Vector::dots(): Vectors do not have compatible Maps:" << std::endl
00119         << "this->getMap(): " << std::endl << *this->getMap()
00120         << "a.getMap(): " << std::endl << *a.getMap() << std::endl);
00121 #else
00122     TEUCHOS_TEST_FOR_EXCEPTION( this->getLocalLength() != a.getLocalLength(), std::runtime_error,
00123         "Tpetra::Vector::dots(): Vectors do not have the same local length.");
00124 #endif
00125     Scalar gbldot;
00126     gbldot = MVT::Dot(this->lclMV_,a.lclMV_);
00127     if (this->isDistributed()) {
00128       Scalar lcldot = gbldot;
00129       Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_SUM,lcldot,outArg(gbldot));
00130     }
00131     return gbldot;
00132   }
00133 
00134   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00135   Scalar Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::meanValue() const {
00136     using Teuchos::as;
00137     using Teuchos::outArg;
00138     typedef Teuchos::ScalarTraits<Scalar> SCT;
00139 
00140     Scalar sum = MVT::Sum(this->lclMV_);
00141     if (this->isDistributed()) {
00142       Scalar lsum = sum;
00143       Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_SUM,lsum,outArg(sum));
00144     }
00145     // mfh 12 Apr 2012: Don't take out the cast from the ordinal type
00146     // to the magnitude type, since operator/ (std::complex<T>, int)
00147     // isn't necessarily defined.
00148     return sum / as<typename SCT::magnitudeType> (this->getGlobalLength ());
00149   }
00150 
00151   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00152   typename Teuchos::ScalarTraits<Scalar>::magnitudeType Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::norm1() const {
00153     using Teuchos::outArg;
00154     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Mag;
00155     Mag norm = MVT::Norm1(this->lclMV_);
00156     if (this->isDistributed()) {
00157       Mag lnorm = norm;
00158       Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_SUM,lnorm,outArg(norm));
00159     }
00160     return norm;
00161   }
00162 
00163   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00164   typename Teuchos::ScalarTraits<Scalar>::magnitudeType Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::norm2() const {
00165     using Teuchos::ScalarTraits;
00166     using Teuchos::outArg;
00167     typedef typename ScalarTraits<Scalar>::magnitudeType Mag;
00168     Mag norm = MVT::Norm2Squared(this->lclMV_);
00169     if (this->isDistributed()) {
00170       Mag lnorm = norm;
00171       Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_SUM,lnorm,outArg(norm));
00172     }
00173     return ScalarTraits<Mag>::squareroot(norm);
00174   }
00175 
00176   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00177   typename Teuchos::ScalarTraits<Scalar>::magnitudeType Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::normInf() const {
00178     using Teuchos::outArg;
00179     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Mag;
00180     Mag norm = MVT::NormInf(this->lclMV_);
00181     if (this->isDistributed()) {
00182       Mag lnorm = norm;
00183       Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_MAX,lnorm,outArg(norm));
00184     }
00185     return norm;
00186   }
00187 
00188   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00189   typename Teuchos::ScalarTraits<Scalar>::magnitudeType Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::normWeighted(const Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &weights) const {
00190     using Teuchos::ScalarTraits;
00191     using Teuchos::outArg;
00192     typedef typename ScalarTraits<Scalar>::magnitudeType Mag;
00193 #ifdef HAVE_TPETRA_DEBUG
00194     TEUCHOS_TEST_FOR_EXCEPTION( !this->getMap()->isCompatible(*weights.getMap()), std::runtime_error,
00195         "Tpetra::Vector::normWeighted(): Vectors do not have compatible Maps:" << std::endl
00196         << "this->getMap(): " << std::endl << *this->getMap()
00197         << "weights.getMap(): " << std::endl << *weights.getMap() << std::endl);
00198 #else
00199     TEUCHOS_TEST_FOR_EXCEPTION( this->getLocalLength() != weights.getLocalLength(), std::runtime_error,
00200         "Tpetra::Vector::normWeighted(): Vectors do not have the same local length.");
00201 #endif
00202     Mag norm = MVT::WeightedNorm(this->lclMV_,weights.lclMV_);
00203     if (this->isDistributed()) {
00204       Mag lnorm = norm;
00205       Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_SUM,lnorm,outArg(norm));
00206     }
00207     return ScalarTraits<Mag>::squareroot(norm / this->getGlobalLength());
00208   }
00209 
00210 
00211   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00212   Teuchos::RCP<const Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
00213   Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
00214   offsetView (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& subMap,
00215               size_t offset) const
00216   {
00217     using Teuchos::RCP;
00218     using Teuchos::rcp;
00219     typedef Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> V;
00220 
00221     const size_t newNumRows = subMap->getNodeNumElements ();
00222     const bool tooManyElts = newNumRows + offset > this->lclMV_.getOrigNumRows ();
00223     if (tooManyElts) {
00224       const int myRank = this->getMap ()->getComm ()->getRank ();
00225       TEUCHOS_TEST_FOR_EXCEPTION(
00226         newNumRows + offset > MVT::getNumRows (this->lclMV_),
00227         std::runtime_error,
00228         "Tpetra::Vector::offsetView: Invalid input Map.  Input Map owns "
00229         << subMap->getNodeNumElements () << " elements on process " << myRank
00230         << ".  offset = " << offset << ".  Yet, the Vector contains only "
00231         << this->lclMV_.getOrigNumRows () << " on this process.");
00232     }
00233 
00234     KokkosClassic::MultiVector<Scalar, Node> newLocalMV =
00235       this->lclMV_.offsetView (newNumRows, this->lclMV_.getNumCols (), offset, 0);
00236     return rcp (new V (subMap, newLocalMV.getValuesNonConst (), COMPUTE_VIEW_CONSTRUCTOR));
00237   }
00238 
00239 
00240   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00241   Teuchos::RCP<Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
00242   Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
00243   offsetViewNonConst (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& subMap,
00244                       size_t offset)
00245   {
00246     using Teuchos::RCP;
00247     using Teuchos::rcp;
00248     typedef Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> V;
00249 
00250     const size_t newNumRows = subMap->getNodeNumElements ();
00251     const bool tooManyElts = newNumRows + offset > this->lclMV_.getOrigNumRows ();
00252     if (tooManyElts) {
00253       const int myRank = this->getMap ()->getComm ()->getRank ();
00254       TEUCHOS_TEST_FOR_EXCEPTION(
00255         newNumRows + offset > MVT::getNumRows (this->lclMV_),
00256         std::runtime_error,
00257         "Tpetra::Vector::offsetViewNonconst: Invalid input Map.  Input Map owns "
00258         << subMap->getNodeNumElements () << " elements on process " << myRank
00259         << ".  offset = " << offset << ".  Yet, the MultiVector contains only "
00260         << this->lclMV_.getOrigNumRows () << " on this process.");
00261     }
00262 
00263     KokkosClassic::MultiVector<Scalar, Node> newLocalMV =
00264       this->lclMV_.offsetViewNonConst (newNumRows, this->lclMV_.getNumCols (), offset, 0);
00265     return rcp (new V (subMap, newLocalMV.getValuesNonConst (), COMPUTE_VIEW_CONSTRUCTOR));
00266   }
00267 
00268 
00269   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00270   std::string
00271   Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::description() const
00272   {
00273     using Teuchos::TypeNameTraits;
00274 
00275     std::ostringstream oss;
00276     oss << "\"Tpetra::Vector\": {"
00277         << "Template parameters: {"
00278         << "Scalar: " << TypeNameTraits<Scalar>::name ()
00279         << ", LocalOrdinal: " << TypeNameTraits<LocalOrdinal>::name ()
00280         << ", GlobalOrdinal: " << TypeNameTraits<GlobalOrdinal>::name ()
00281         << ", Node: " << TypeNameTraits<Node>::name ()
00282         << "}";
00283     if (this->getObjectLabel () != "") {
00284       oss << ", Label: \"" << this->getObjectLabel () << "\", ";
00285     }
00286     oss << "Global length: " << this->getGlobalLength ()
00287         << "}";
00288     return oss.str ();
00289   }
00290 
00291   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00292   void Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
00293   describe (Teuchos::FancyOStream& out,
00294             const Teuchos::EVerbosityLevel verbLevel) const
00295   {
00296     using std::endl;
00297     using std::setw;
00298     using Teuchos::Comm;
00299     using Teuchos::RCP;
00300     using Teuchos::TypeNameTraits;
00301     using Teuchos::VERB_DEFAULT;
00302     using Teuchos::VERB_NONE;
00303     using Teuchos::VERB_LOW;
00304     using Teuchos::VERB_MEDIUM;
00305     using Teuchos::VERB_HIGH;
00306     using Teuchos::VERB_EXTREME;
00307 
00308     const Teuchos::EVerbosityLevel vl =
00309       (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
00310 
00311     const Map<LocalOrdinal, GlobalOrdinal, Node>& map = * (this->getMap ());
00312     RCP<const Comm<int> > comm = map.getComm ();
00313     const int myImageID = comm->getRank ();
00314     const int numImages = comm->getSize ();
00315     Teuchos::OSTab tab0 (out);
00316 
00317     if (vl != VERB_NONE) {
00318       if (myImageID == 0) {
00319         out << "\"Tpetra::Vector\":" << endl;
00320       }
00321       Teuchos::OSTab tab1 (out);// applies to all processes
00322       if (myImageID == 0) {
00323         out << "Template parameters:";
00324         {
00325           Teuchos::OSTab tab2 (out);
00326           out << "Scalar: " << TypeNameTraits<Scalar>::name () << endl
00327               << "LocalOrdinal: " << TypeNameTraits<LocalOrdinal>::name () << endl
00328               << "GlobalOrdinal: " << TypeNameTraits<GlobalOrdinal>::name () << endl
00329               << "Node: " << TypeNameTraits<Node>::name () << endl;
00330         }
00331         out << endl;
00332         if (this->getObjectLabel () != "") {
00333           out << "Label: \"" << this->getObjectLabel () << "\"" << endl;
00334         }
00335         out << "Global length: " << this->getGlobalLength () << endl;
00336       }
00337       for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
00338         if (myImageID == imageCtr) {
00339           if (vl != VERB_LOW) {
00340             // VERB_MEDIUM and higher prints getLocalLength()
00341             out << "Process: " << myImageID << endl;
00342             Teuchos::OSTab tab2 (out);
00343             out << "Local length: " << this->getLocalLength () << endl;
00344 
00345             if (vl == VERB_EXTREME && this->getLocalLength () > 0) {
00346               // VERB_EXTREME prints values
00347               out << "Global indices and values:" << endl;
00348               Teuchos::OSTab tab3 (out);
00349               RCP<Node> node = this->lclMV_.getNode ();
00350               ArrayRCP<const Scalar> myview =
00351                 node->template viewBuffer<Scalar> (this->getLocalLength (),
00352                                                    MVT::getValues (this->lclMV_));
00353               for (size_t i = 0; i < this->getLocalLength (); ++i) {
00354                 out << map.getGlobalElement (i) << ": " << myview[i] << endl;
00355               }
00356             }
00357           }
00358           std::flush (out); // give output time to complete
00359         }
00360         comm->barrier (); // give output time to complete
00361         comm->barrier ();
00362         comm->barrier ();
00363       }
00364     }
00365   }
00366 
00367   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00368   Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node >
00369     createCopy( const Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node >& src) {
00370     return src;
00371   }
00372 } // namespace Tpetra
00373 
00374 //
00375 // Explicit instantiation macro
00376 //
00377 // Must be expanded from within the Tpetra namespace!
00378 //
00379 
00380 #if defined(TPETRA_HAVE_KOKKOS_REFACTOR)
00381 #include "Tpetra_KokkosRefactor_Vector_def.hpp"
00382 #endif
00383 
00384 
00385 #define TPETRA_VECTOR_INSTANT(SCALAR,LO,GO,NODE) \
00386   \
00387   template class Vector< SCALAR , LO , GO , NODE >; \
00388   template Vector< SCALAR , LO , GO , NODE > createCopy( const Vector< SCALAR , LO , GO , NODE >& src); \
00389 
00390 
00391 #endif // TPETRA_VECTOR_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines