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

Generated on Tue Jul 13 09:39:07 2010 for Tpetra Matrix/Vector Services by  doxygen 1.4.7