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