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 <Teuchos_ScalarTraits.hpp>
00035 #include <Teuchos_OrdinalTraits.hpp>
00036 #include <Teuchos_TypeNameTraits.hpp>
00037 #include <Teuchos_Tuple.hpp>
00038 
00039 #include "Tpetra_MultiVector.hpp"
00040 
00041 #ifdef DOXYGEN_USE_ONLY
00042   #include "Tpetra_Vector_decl.hpp"
00043 #endif
00044 
00045 namespace Tpetra {
00046 
00047   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00048   Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Vector(const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map, bool zeroOut) 
00049     : MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map,1,zeroOut) {
00050   }
00051 
00052   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00053   Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Vector(const Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &source)
00054   : MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(source) {
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, const Teuchos::ArrayView<const Scalar> &values)
00059   : MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map,values,values.size(),1) {
00060   }
00061 
00062   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00063   Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Vector(const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map, Teuchos::ArrayRCP<Scalar> values)
00064     : MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map,values,map->getNodeNumElements(),1) {
00065   }
00066 
00067   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00068   Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::~Vector() {}
00069 
00070   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00071   void Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::replaceGlobalValue(GlobalOrdinal globalRow, const Scalar &value) {
00072     this->MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::replaceGlobalValue(globalRow,0,value);
00073   }
00074 
00075   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00076   void Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::sumIntoGlobalValue(GlobalOrdinal globalRow, const Scalar &value) {
00077     this->MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::sumIntoGlobalValue(globalRow,0,value);
00078   }
00079 
00080   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00081   void Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::replaceLocalValue(LocalOrdinal myRow, const Scalar &value) {
00082     this->MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::replaceLocalValue(myRow,0,value);
00083   }
00084 
00085   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00086   void Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::sumIntoLocalValue(LocalOrdinal myRow, const Scalar &value) {
00087     this->MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::sumIntoLocalValue(myRow,0,value);
00088   }
00089 
00090   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00091   void Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::get1dCopy(Teuchos::ArrayView<Scalar> A) const {
00092     size_t lda = this->getLocalLength();
00093     this->get1dCopy(A,lda);
00094   }
00095 
00096   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00097   Scalar Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::dot(const Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &a) const {
00098     using Teuchos::outArg;
00099 #ifdef HAVE_TPETRA_DEBUG
00100     TEST_FOR_EXCEPTION( !this->getMap()->isCompatible(*a.getMap()), std::runtime_error,
00101         "Tpetra::Vector::dots(): Vectors do not have compatible Maps:" << std::endl
00102         << "this->getMap(): " << std::endl << *this->getMap() 
00103         << "a.getMap(): " << std::endl << *a.getMap() << std::endl);
00104 #else
00105     TEST_FOR_EXCEPTION( this->getLocalLength() != a.getLocalLength(), std::runtime_error,
00106         "Tpetra::Vector::dots(): Vectors do not have the same local length.");
00107 #endif
00108     Scalar gbldot;
00109     gbldot = MVT::Dot(this->lclMV_,a.lclMV_);
00110     if (this->isDistributed()) {
00111       Scalar lcldot = gbldot;
00112       Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_SUM,lcldot,outArg(gbldot));
00113     }
00114     return gbldot;
00115   }
00116 
00117   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00118   Scalar Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::meanValue() const {
00119     using Teuchos::outArg;
00120     typedef Teuchos::ScalarTraits<Scalar> SCT;
00121     Scalar sum = MVT::Sum(this->lclMV_);
00122     if (this->isDistributed()) {
00123       Scalar lsum = sum;
00124       Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_SUM,lsum,outArg(sum));
00125     }
00126     return sum / Teuchos::as<Scalar>(this->getGlobalLength());
00127   }
00128 
00129   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00130   typename Teuchos::ScalarTraits<Scalar>::magnitudeType Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::norm1() const {
00131     using Teuchos::outArg;
00132     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Mag;
00133     Mag norm = MVT::Norm1(this->lclMV_);
00134     if (this->isDistributed()) {
00135       Mag lnorm = norm;
00136       Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_SUM,lnorm,outArg(norm));
00137     }
00138     return norm;
00139   }
00140 
00141   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00142   typename Teuchos::ScalarTraits<Scalar>::magnitudeType Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::norm2() const {
00143     using Teuchos::ScalarTraits;
00144     using Teuchos::outArg;
00145     typedef typename ScalarTraits<Scalar>::magnitudeType Mag;
00146     Mag norm = MVT::Norm2Squared(this->lclMV_);
00147     if (this->isDistributed()) {
00148       Mag lnorm = norm;
00149       Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_SUM,lnorm,outArg(norm));
00150     }
00151     return ScalarTraits<Mag>::squareroot(norm);
00152   }
00153 
00154   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00155   typename Teuchos::ScalarTraits<Scalar>::magnitudeType Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::normInf() const {
00156     using Teuchos::outArg;
00157     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Mag;
00158     Mag norm = MVT::NormInf(this->lclMV_);
00159     if (this->isDistributed()) {
00160       Mag lnorm = norm;
00161       Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_MAX,lnorm,outArg(norm));
00162     }
00163     return norm;
00164   }
00165 
00166   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00167   typename Teuchos::ScalarTraits<Scalar>::magnitudeType Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::normWeighted(const Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &weights) const {
00168     using Teuchos::ScalarTraits;
00169     using Teuchos::outArg;
00170     typedef typename ScalarTraits<Scalar>::magnitudeType Mag;
00171 #ifdef HAVE_TPETRA_DEBUG
00172     TEST_FOR_EXCEPTION( !this->getMap()->isCompatible(*weights.getMap()), std::runtime_error,
00173         "Tpetra::Vector::normWeighted(): Vectors do not have compatible Maps:" << std::endl
00174         << "this->getMap(): " << std::endl << *this->getMap() 
00175         << "weights.getMap(): " << std::endl << *weights.getMap() << std::endl);
00176 #else
00177     TEST_FOR_EXCEPTION( this->getLocalLength() != weights.getLocalLength(), std::runtime_error,
00178         "Tpetra::Vector::normWeighted(): Vectors do not have the same local length.");
00179 #endif
00180     Mag norm = MVT::WeightedNorm(this->lclMV_,weights.lclMV_);
00181     if (this->isDistributed()) {
00182       Mag lnorm = norm;
00183       Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_SUM,lnorm,outArg(norm));
00184     }
00185     return ScalarTraits<Mag>::squareroot(norm / Teuchos::as<Mag>(this->getGlobalLength()));
00186   }
00187 
00188   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00189   std::string Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::description() const {
00190     std::ostringstream oss;
00191     oss << Teuchos::Describable::description();
00192     oss << "{length="<<this->getGlobalLength()
00193         << "}";
00194     return oss.str();
00195   }
00196 
00197   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00198   void Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const {
00199     using std::endl;
00200     using std::setw;
00201     using Teuchos::VERB_DEFAULT;
00202     using Teuchos::VERB_NONE;
00203     using Teuchos::VERB_LOW;
00204     using Teuchos::VERB_MEDIUM;
00205     using Teuchos::VERB_HIGH;
00206     using Teuchos::VERB_EXTREME;
00207     Teuchos::EVerbosityLevel vl = verbLevel;
00208     if (vl == VERB_DEFAULT) vl = VERB_LOW;
00209     Teuchos::RCP<const Teuchos::Comm<int> > comm = this->getMap()->getComm();
00210     const int myImageID = comm->getRank(),
00211               numImages = comm->getSize();
00212     size_t width = 1;
00213     for (size_t dec=10; dec<this->getGlobalLength(); dec *= 10) {
00214       ++width;
00215     }
00216     Teuchos::OSTab tab(out);
00217     if (vl != VERB_NONE) {
00218       // VERB_LOW and higher prints description()
00219       if (myImageID == 0) out << this->description() << std::endl; 
00220       for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
00221         if (myImageID == imageCtr) {
00222           if (vl != VERB_LOW) {
00223             // VERB_MEDIUM and higher prints getLocalLength()
00224             out << "node " << setw(width) << myImageID << ": local length=" << this->getLocalLength() << endl;
00225             if (vl != VERB_MEDIUM) {
00226               // VERB_HIGH and higher prints isConstantStride() and stride()
00227               if (vl == VERB_EXTREME && this->getLocalLength() > 0) {
00228                 Teuchos::RCP<Node> node = this->lclMV_.getNode();
00229                 KOKKOS_NODE_TRACE("Vector::describe()")
00230                 Teuchos::ArrayRCP<const Scalar> myview = node->template viewBuffer<Scalar>(
00231                                                                this->getLocalLength(), 
00232                                                                MVT::getValues(this->lclMV_) );
00233                 // VERB_EXTREME prints values
00234                 for (size_t i=0; i<this->getLocalLength(); ++i) {
00235                   out << setw(width) << this->getMap()->getGlobalElement(i) 
00236                       << ": "
00237                       << myview[i] << endl;
00238                 }
00239                 myview = Teuchos::null;
00240               }
00241             }
00242             else {
00243               out << endl;
00244             }
00245           }
00246         }
00247       }
00248     }
00249   }
00250 
00251 } // namespace Tpetra
00252 
00253 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00254 Teuchos::RCP< Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
00255 Tpetra::createVector(const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > &map) {
00256   const bool DO_INIT_TO_ZERO = true;
00257   return Teuchos::rcp( new Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map,DO_INIT_TO_ZERO) );
00258 }
00259 
00260 //
00261 // Explicit instantiation macro
00262 //
00263 // Must be expanded from within the Tpetra namespace!
00264 //
00265 
00266 #define TPETRA_VECTOR_INSTANT(SCALAR,LO,GO,NODE) \
00267   \
00268   template class Vector< SCALAR , LO , GO , NODE >; \
00269   \
00270   template Teuchos::RCP< Vector<SCALAR,LO,GO,NODE> > \
00271   createVector<SCALAR,LO,GO,NODE>(const Teuchos::RCP< const Tpetra::Map<LO,GO,NODE> > &map);
00272 
00273 
00274 #endif // TPETRA_VECTOR_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:21:41 2011 for Tpetra Matrix/Vector Services by  doxygen 1.6.3