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 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   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00211   std::string Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::description() const {
00212     std::ostringstream oss;
00213     oss << Teuchos::Describable::description();
00214     oss << "{length="<<this->getGlobalLength()
00215         << "}";
00216     return oss.str();
00217   }
00218 
00219   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00220   void Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const {
00221     using std::endl;
00222     using std::setw;
00223     using Teuchos::VERB_DEFAULT;
00224     using Teuchos::VERB_NONE;
00225     using Teuchos::VERB_LOW;
00226     using Teuchos::VERB_MEDIUM;
00227     using Teuchos::VERB_HIGH;
00228     using Teuchos::VERB_EXTREME;
00229     Teuchos::EVerbosityLevel vl = verbLevel;
00230     if (vl == VERB_DEFAULT) vl = VERB_LOW;
00231     RCP<const Teuchos::Comm<int> > comm = this->getMap()->getComm();
00232     const int myImageID = comm->getRank(),
00233               numImages = comm->getSize();
00234     size_t width = 1;
00235     for (size_t dec=10; dec<this->getGlobalLength(); dec *= 10) {
00236       ++width;
00237     }
00238     Teuchos::OSTab tab(out);
00239     if (vl != VERB_NONE) {
00240       // VERB_LOW and higher prints description()
00241       if (myImageID == 0) out << this->description() << std::endl; 
00242       for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
00243         if (myImageID == imageCtr) {
00244           if (vl != VERB_LOW) {
00245             // VERB_MEDIUM and higher prints getLocalLength()
00246             out << "node " << setw(width) << myImageID << ": local length=" << this->getLocalLength() << endl;
00247             if (vl != VERB_MEDIUM) {
00248               // VERB_HIGH and higher prints isConstantStride() and stride()
00249               if (vl == VERB_EXTREME && this->getLocalLength() > 0) {
00250                 RCP<Node> node = this->lclMV_.getNode();
00251                 KOKKOS_NODE_TRACE("Vector::describe()")
00252                 ArrayRCP<const Scalar> myview = node->template viewBuffer<Scalar>(
00253                                                                this->getLocalLength(), 
00254                                                                MVT::getValues(this->lclMV_) );
00255                 // VERB_EXTREME prints values
00256                 for (size_t i=0; i<this->getLocalLength(); ++i) {
00257                   out << setw(width) << this->getMap()->getGlobalElement(i) 
00258                       << ": "
00259                       << myview[i] << endl;
00260                 }
00261                 myview = Teuchos::null;
00262               }
00263             }
00264             else {
00265               out << endl;
00266             }
00267           }
00268         }
00269       }
00270     }
00271   }
00272 
00273 } // namespace Tpetra
00274 
00275 //
00276 // Explicit instantiation macro
00277 //
00278 // Must be expanded from within the Tpetra namespace!
00279 //
00280 
00281 #define TPETRA_VECTOR_INSTANT(SCALAR,LO,GO,NODE) \
00282   \
00283   template class Vector< SCALAR , LO , GO , NODE >; \
00284 
00285 
00286 #endif // TPETRA_VECTOR_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines