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 terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef TPETRA_VECTOR_DEF_HPP
00030 #define TPETRA_VECTOR_DEF_HPP
00031 
00032 #include <Kokkos_NodeTrace.hpp>
00033 
00034 #include "Tpetra_MultiVector.hpp"
00035 
00036 #ifdef DOXYGEN_USE_ONLY
00037   #include "Tpetra_Vector_decl.hpp"
00038 #endif
00039 
00040 namespace Tpetra {
00041 
00042   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00043   Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Vector(const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map, bool zeroOut) 
00044     : MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map,1,zeroOut) {
00045   }
00046 
00047   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00048   Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Vector(const Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &source)
00049   : MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(source) {
00050   }
00051 
00052   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00053   Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Vector(const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map, const Teuchos::ArrayView<const Scalar> &values)
00054   : MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map,values,values.size(),1) {
00055   }
00056 
00057   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00058   Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Vector(const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map, Teuchos::ArrayRCP<Scalar> values)
00059     : MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map,values,map->getNodeNumElements(),1) {
00060   }
00061 
00062   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00063   Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::~Vector() {}
00064 
00065   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00066   void Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::replaceGlobalValue(GlobalOrdinal globalRow, const Scalar &value) {
00067     this->MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::replaceGlobalValue(globalRow,0,value);
00068   }
00069 
00070   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00071   void Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::sumIntoGlobalValue(GlobalOrdinal globalRow, const Scalar &value) {
00072     this->MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::sumIntoGlobalValue(globalRow,0,value);
00073   }
00074 
00075   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00076   void Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::replaceLocalValue(LocalOrdinal myRow, const Scalar &value) {
00077     this->MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::replaceLocalValue(myRow,0,value);
00078   }
00079 
00080   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00081   void Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::sumIntoLocalValue(LocalOrdinal myRow, const Scalar &value) {
00082     this->MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::sumIntoLocalValue(myRow,0,value);
00083   }
00084 
00085   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00086   void Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::get1dCopy(Teuchos::ArrayView<Scalar> A) const {
00087     size_t lda = this->getLocalLength();
00088     this->get1dCopy(A,lda);
00089   }
00090 
00091   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00092   Scalar Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::dot(const Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &a) const {
00093     using Teuchos::outArg;
00094 #ifdef HAVE_TPETRA_DEBUG
00095     TEST_FOR_EXCEPTION( !this->getMap()->isCompatible(*a.getMap()), std::runtime_error,
00096         "Tpetra::Vector::dots(): Vectors do not have compatible Maps:" << std::endl
00097         << "this->getMap(): " << std::endl << *this->getMap() 
00098         << "a.getMap(): " << std::endl << *a.getMap() << std::endl);
00099 #else
00100     TEST_FOR_EXCEPTION( this->getLocalLength() != a.getLocalLength(), std::runtime_error,
00101         "Tpetra::Vector::dots(): Vectors do not have the same local length.");
00102 #endif
00103     Scalar gbldot;
00104     gbldot = MVT::Dot(this->lclMV_,a.lclMV_);
00105     if (this->isDistributed()) {
00106       Scalar lcldot = gbldot;
00107       Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_SUM,lcldot,outArg(gbldot));
00108     }
00109     return gbldot;
00110   }
00111 
00112   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00113   Scalar Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::meanValue() const {
00114     using Teuchos::outArg;
00115     typedef Teuchos::ScalarTraits<Scalar> SCT;
00116     Scalar sum = MVT::Sum(this->lclMV_);
00117     if (this->isDistributed()) {
00118       Scalar lsum = sum;
00119       Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_SUM,lsum,outArg(sum));
00120     }
00121     return sum / Teuchos::as<Scalar>(this->getGlobalLength());
00122   }
00123 
00124   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00125   typename Teuchos::ScalarTraits<Scalar>::magnitudeType Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::norm1() const {
00126     using Teuchos::outArg;
00127     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Mag;
00128     Mag norm = MVT::Norm1(this->lclMV_);
00129     if (this->isDistributed()) {
00130       Mag lnorm = norm;
00131       Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_SUM,lnorm,outArg(norm));
00132     }
00133     return norm;
00134   }
00135 
00136   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00137   typename Teuchos::ScalarTraits<Scalar>::magnitudeType Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::norm2() const {
00138     using Teuchos::ScalarTraits;
00139     using Teuchos::outArg;
00140     typedef typename ScalarTraits<Scalar>::magnitudeType Mag;
00141     Mag norm = MVT::Norm2Squared(this->lclMV_);
00142     if (this->isDistributed()) {
00143       Mag lnorm = norm;
00144       Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_SUM,lnorm,outArg(norm));
00145     }
00146     return ScalarTraits<Mag>::squareroot(norm);
00147   }
00148 
00149   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00150   typename Teuchos::ScalarTraits<Scalar>::magnitudeType Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::normInf() const {
00151     using Teuchos::outArg;
00152     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Mag;
00153     Mag norm = MVT::NormInf(this->lclMV_);
00154     if (this->isDistributed()) {
00155       Mag lnorm = norm;
00156       Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_MAX,lnorm,outArg(norm));
00157     }
00158     return norm;
00159   }
00160 
00161   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00162   typename Teuchos::ScalarTraits<Scalar>::magnitudeType Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::normWeighted(const Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &weights) const {
00163     using Teuchos::ScalarTraits;
00164     using Teuchos::outArg;
00165     typedef typename ScalarTraits<Scalar>::magnitudeType Mag;
00166 #ifdef HAVE_TPETRA_DEBUG
00167     TEST_FOR_EXCEPTION( !this->getMap()->isCompatible(*weights.getMap()), std::runtime_error,
00168         "Tpetra::Vector::normWeighted(): Vectors do not have compatible Maps:" << std::endl
00169         << "this->getMap(): " << std::endl << *this->getMap() 
00170         << "weights.getMap(): " << std::endl << *weights.getMap() << std::endl);
00171 #else
00172     TEST_FOR_EXCEPTION( this->getLocalLength() != weights.getLocalLength(), std::runtime_error,
00173         "Tpetra::Vector::normWeighted(): Vectors do not have the same local length.");
00174 #endif
00175     Mag norm = MVT::WeightedNorm(this->lclMV_,weights.lclMV_);
00176     if (this->isDistributed()) {
00177       Mag lnorm = norm;
00178       Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_SUM,lnorm,outArg(norm));
00179     }
00180     return ScalarTraits<Mag>::squareroot(norm / Teuchos::as<Mag>(this->getGlobalLength()));
00181   }
00182 
00183   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00184   std::string Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::description() const {
00185     std::ostringstream oss;
00186     oss << Teuchos::Describable::description();
00187     oss << "{length="<<this->getGlobalLength()
00188         << "}";
00189     return oss.str();
00190   }
00191 
00192   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00193   void Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const {
00194     using std::endl;
00195     using std::setw;
00196     using Teuchos::VERB_DEFAULT;
00197     using Teuchos::VERB_NONE;
00198     using Teuchos::VERB_LOW;
00199     using Teuchos::VERB_MEDIUM;
00200     using Teuchos::VERB_HIGH;
00201     using Teuchos::VERB_EXTREME;
00202     Teuchos::EVerbosityLevel vl = verbLevel;
00203     if (vl == VERB_DEFAULT) vl = VERB_LOW;
00204     Teuchos::RCP<const Teuchos::Comm<int> > comm = this->getMap()->getComm();
00205     const int myImageID = comm->getRank(),
00206               numImages = comm->getSize();
00207     size_t width = 1;
00208     for (size_t dec=10; dec<this->getGlobalLength(); dec *= 10) {
00209       ++width;
00210     }
00211     Teuchos::OSTab tab(out);
00212     if (vl != VERB_NONE) {
00213       // VERB_LOW and higher prints description()
00214       if (myImageID == 0) out << this->description() << std::endl; 
00215       for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
00216         if (myImageID == imageCtr) {
00217           if (vl != VERB_LOW) {
00218             // VERB_MEDIUM and higher prints getLocalLength()
00219             out << "node " << setw(width) << myImageID << ": local length=" << this->getLocalLength() << endl;
00220             if (vl != VERB_MEDIUM) {
00221               // VERB_HIGH and higher prints isConstantStride() and stride()
00222               if (vl == VERB_EXTREME && this->getLocalLength() > 0) {
00223                 Teuchos::RCP<Node> node = this->lclMV_.getNode();
00224                 KOKKOS_NODE_TRACE("Vector::describe()")
00225                 Teuchos::ArrayRCP<const Scalar> myview = node->template viewBuffer<Scalar>(
00226                                                                this->getLocalLength(), 
00227                                                                MVT::getValues(this->lclMV_) );
00228                 // VERB_EXTREME prints values
00229                 for (size_t i=0; i<this->getLocalLength(); ++i) {
00230                   out << setw(width) << this->getMap()->getGlobalElement(i) 
00231                       << ": "
00232                       << myview[i] << endl;
00233                 }
00234                 myview = Teuchos::null;
00235               }
00236             }
00237             else {
00238               out << endl;
00239             }
00240           }
00241         }
00242       }
00243     }
00244   }
00245 
00246 } // namespace Tpetra
00247 
00248 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00249 Teuchos::RCP< Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
00250 Tpetra::createVector(const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > &map) {
00251   const bool DO_INIT_TO_ZERO = true;
00252   return Teuchos::rcp( new Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map,DO_INIT_TO_ZERO) );
00253 }
00254 
00255 //
00256 // Explicit instantiation macro
00257 //
00258 // Must be expanded from within the Tpetra namespace!
00259 //
00260 
00261 #define TPETRA_VECTOR_INSTANT(SCALAR,LO,GO,NODE) \
00262   \
00263   template class Vector< SCALAR , LO , GO , NODE >; \
00264   \
00265   template Teuchos::RCP< Vector<SCALAR,LO,GO,NODE> > \
00266   createVector<SCALAR,LO,GO,NODE>(const Teuchos::RCP< const Tpetra::Map<LO,GO,NODE> > &map);
00267 
00268 
00269 #endif // TPETRA_VECTOR_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines