Tpetra Matrix/Vector Services: Tpetra

00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //          Tpetra: Templated Linear Algebra Services Package
00005 //                 Copyright (2004) 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_HPP
00030 #define TPETRA_VECTOR_HPP
00031 
00032 #include <Teuchos_ScalarTraits.hpp>
00033 #include <Teuchos_OrdinalTraits.hpp>
00034 #include <Teuchos_Object.hpp>
00035 #include <Teuchos_BLAS.hpp>
00036 
00037 #include "Tpetra_VectorDecl.hpp"
00038 #include "Tpetra_VectorData.hpp"
00039 
00040 namespace Tpetra {
00041 
00042   template<typename Ordinal,typename Scalar>
00043   Vector<Ordinal,Scalar>::Vector(const Map<Ordinal> &map) 
00044     : DistObject<Ordinal,Scalar>(map, map.getPlatform()->createComm(), "Tpetra::Vector")
00045     , VectorData_()
00046   {
00047     const Scalar ZERO = Teuchos::ScalarTraits<Scalar>::zero();
00048     VectorData_ = Teuchos::rcp(new VectorData<Ordinal,Scalar>(map, ZERO));
00049   }
00050 
00051   template<typename Ordinal,typename Scalar>
00052   Vector<Ordinal,Scalar>::Vector(const Teuchos::ArrayView<const Scalar> &values, const Map<Ordinal> &map)
00053   : DistObject<Ordinal,Scalar>(map,map.getPlatform()->createComm(), "Tpetra::Vector")
00054   , VectorData_()
00055   {
00056     const Scalar ZERO = Teuchos::ScalarTraits<Scalar>::zero();
00057     const Ordinal OZERO = Teuchos::OrdinalTraits<Ordinal>::zero();
00058     const Ordinal length = map.getNumMyEntries();
00059 
00060     VectorData_ = Teuchos::rcp(new VectorData<Ordinal,Scalar>(map, ZERO));
00061 
00062     TEST_FOR_EXCEPTION(values.size() != length, std::invalid_argument,
00063       "Tpetra::Vector::constructor(values,map): number of values does not match size of map.");
00064     std::copy(values.begin(),values.end(), VectorData_->values_().begin());
00065   }
00066 
00067   template<typename Ordinal, typename Scalar>
00068   Vector<Ordinal,Scalar>::Vector(const Vector<Ordinal,Scalar> &source)
00069   : DistObject<Ordinal,Scalar>(source)
00070   , VectorData_(source.VectorData_)
00071   {}
00072 
00073   template<typename Ordinal, typename Scalar>
00074   Vector<Ordinal,Scalar>::~Vector() {}
00075 
00076   template<typename Ordinal, typename Scalar>
00077   void Vector<Ordinal,Scalar>::submitEntries(const Teuchos::ArrayView<const Ordinal> &indices,
00078                                              const Teuchos::ArrayView<const Scalar>  &values)
00079   {
00080     TEST_FOR_EXCEPTION(indices.size() != values.size(), std::invalid_argument,
00081         "Tpetra::Vector::submitEntries(indices,values): finish.");
00082     const Ordinal ordinalZero = Teuchos::OrdinalTraits<Ordinal>::zero();
00083     typename Teuchos::ArrayView<const Ordinal>::const_iterator ind = indices.begin();
00084     typename Teuchos::ArrayView<const Scalar>::const_iterator  val = values.begin();
00085     for (; ind != indices.end(); ++ind, ++val) {
00086       // if TEUCHOS_DEBUG, this is bounds checked
00087       VectorData_->values_[*ind] += *val;
00088     }
00089   }
00090 
00091   template<typename Ordinal, typename Scalar>
00092   void Vector<Ordinal,Scalar>::setAllToScalar(const Scalar &value) {
00093     TEST_FOR_EXCEPT(true); // finish
00094   }
00095 
00096   template<typename Ordinal, typename Scalar>
00097   void Vector<Ordinal,Scalar>::setAllToRandom() {
00098     TEST_FOR_EXCEPT(true); // finish
00099   }
00100 
00101   template<typename Ordinal, typename Scalar>
00102   Scalar Vector<Ordinal,Scalar>::dotProduct(const Vector<Ordinal,Scalar> &x) const 
00103   {
00104     TEST_FOR_EXCEPTION( !getMap().isCompatible(x.getMap()), 
00105         std::runtime_error, "Tpetra::Vector::dotProduct(): Vectors must have compatible Maps.");
00106     Teuchos::BLAS<Ordinal,Scalar> blas;
00107     // call BLAS routine to calculate local dot product
00108     // use Comm call to sum all local dot products
00109     // update flops counter
00110     TEST_FOR_EXCEPT(true);
00111     return Teuchos::ScalarTraits<Scalar>::zero();
00112   }
00113 
00114   template<typename Ordinal, typename Scalar>
00115   void Vector<Ordinal,Scalar>::absoluteValue(const Vector<Ordinal,Scalar> &x) 
00116   {
00117     TEST_FOR_EXCEPTION( !getMap().isCompatible(x.getMap()), 
00118         std::runtime_error, "Tpetra::Vector::absoluteValue(): Vectors must have compatible Maps.");
00119     // do work
00120     // update flops counter
00121     TEST_FOR_EXCEPT(true);
00122   }
00123 
00124   template<typename Ordinal, typename Scalar>
00125   void Vector<Ordinal,Scalar>::reciprocal(const Vector<Ordinal,Scalar> &x) 
00126   {
00127     TEST_FOR_EXCEPTION( !getMap().isCompatible(x.getMap()), 
00128         std::runtime_error, "Tpetra::Vector::reciprocal(): Vectors must have compatible Maps.");
00129     // do work
00130     // update flops counter
00131     TEST_FOR_EXCEPT(true);
00132   }
00133 
00134   template<typename Ordinal, typename Scalar>
00135   void Vector<Ordinal,Scalar>::scale(const Scalar &scalarThis) 
00136   {
00137     // do work
00138     // update flops counter
00139     TEST_FOR_EXCEPT(true);
00140   }
00141 
00142   template<typename Ordinal, typename Scalar>
00143   void Vector<Ordinal,Scalar>::scale(const Scalar &scalarX, const Vector<Ordinal,Scalar> &x) 
00144   {
00145     TEST_FOR_EXCEPTION( !getMap().isCompatible(x.getMap()), 
00146         std::runtime_error, "Tpetra::Vector::scale(): Vectors must have compatible Maps.");
00147     // do work
00148     // update flops counter
00149     TEST_FOR_EXCEPT(true);
00150   }
00151 
00152   template<typename Ordinal, typename Scalar>
00153   void Vector<Ordinal,Scalar>::update(const Scalar &scalarX, const Vector<Ordinal,Scalar> &x, 
00154                                       const Scalar &scalarThis) 
00155   {
00156     TEST_FOR_EXCEPTION( !getMap().isCompatible(x.getMap()), 
00157         std::runtime_error, "Tpetra::Vector::update(): Vectors must have compatible Maps.");
00158     // do work
00159     // update flops counter
00160     TEST_FOR_EXCEPT(true);
00161   }
00162 
00163 
00164   template<typename Ordinal, typename Scalar>
00165   void Vector<Ordinal,Scalar>::update(const Scalar &scalarX, const Vector<Ordinal,Scalar> &x, 
00166                                       const Scalar &scalarY, const Vector<Ordinal,Scalar> &y, 
00167                                       const Scalar &scalarThis) 
00168   {
00169     TEST_FOR_EXCEPTION( !getMap().isCompatible(x.getMap()), 
00170         std::runtime_error, "Tpetra::Vector::update(x,y): Vectors must have compatible Maps.");
00171     TEST_FOR_EXCEPTION( !getMap().isCompatible(y.getMap()), 
00172         std::runtime_error, "Tpetra::Vector::update(x,y): Vectors must have compatible Maps.");
00173     // do work
00174     // update flops counter
00175     TEST_FOR_EXCEPT(true);
00176   }
00177 
00178 
00179   template<typename Ordinal, typename Scalar>
00180   Scalar Vector<Ordinal,Scalar>::norm1() const {
00181     // do work
00182     // update flops counter
00183     TEST_FOR_EXCEPT(true);
00184   }
00185 
00186 
00187   template<typename Ordinal, typename Scalar>
00188   Scalar Vector<Ordinal,Scalar>::norm2() const {
00189     // do work
00190     // update flops counter
00191     TEST_FOR_EXCEPT(true);
00192   }
00193 
00194 
00195   template<typename Ordinal, typename Scalar>
00196   Scalar Vector<Ordinal,Scalar>::normInf() const {
00197     // do work
00198     // update flops counter
00199     TEST_FOR_EXCEPT(true);
00200   }
00201 
00202 
00203   template<typename Ordinal, typename Scalar>
00204   Scalar Vector<Ordinal,Scalar>::normWeighted(const Vector<Ordinal,Scalar> &weights) const 
00205   {
00206     TEST_FOR_EXCEPTION( !getMap().isCompatible(weights.getMap()), 
00207         std::runtime_error, "Tpetra::Vector::normWeighted(): Vectors must have compatible Maps.");
00208     // do work
00209     // update flops counter
00210     TEST_FOR_EXCEPT(true);
00211   }
00212 
00213 
00214   // finish: move these somewhere?
00215   /*
00216   template<typename Ordinal, typename Scalar>
00217   struct vector_less_mag : public binary_function<Scalar,Scalar,bool> {
00218     bool operator()(Scalar x, Scalar y) { 
00219       return Teuchos::ScalarTraits<Scalar>::magnitude(x) < 
00220         Teuchos::ScalarTraits<Scalar>::magnitude(y); 
00221     }
00222   };
00223 
00224   template<typename Ordinal, typename Scalar>
00225   struct vector_greater_mag : public binary_function<Scalar, Scalar, bool> {
00226     bool operator()(Scalar x, Scalar y) { 
00227       return Teuchos::ScalarTraits<Scalar>::magnitude(x) > 
00228         Teuchos::ScalarTraits<Scalar>::magnitude(y); 
00229     }
00230   };
00231   */
00232 
00233 
00234   template<typename Ordinal, typename Scalar>
00235   Scalar Vector<Ordinal,Scalar>::minValue() const {
00236     TEST_FOR_EXCEPT(true);
00237     return Teuchos::ScalarTraits<Scalar>::zero();
00238   }
00239 
00240 
00241   template<typename Ordinal, typename Scalar>
00242   Scalar Vector<Ordinal,Scalar>::maxValue() const {
00243     TEST_FOR_EXCEPT(true);
00244     return Teuchos::ScalarTraits<Scalar>::zero();
00245   }
00246 
00247 
00248   template<typename Ordinal, typename Scalar>
00249   Scalar Vector<Ordinal,Scalar>::meanValue() const {
00250     TEST_FOR_EXCEPT(true);
00251     // update flops
00252     return Teuchos::ScalarTraits<Scalar>::zero();
00253   }
00254 
00255 
00256   template<typename Ordinal, typename Scalar>
00257   void Vector<Ordinal,Scalar>::elementwiseMultiply(const Scalar &scalarXY, const Vector<Ordinal,Scalar> &x, const Vector<Ordinal,Scalar> &y, 
00258                                                    const Scalar &scalarThis)
00259   {
00260     TEST_FOR_EXCEPTION( !getMap().isCompatible(x.getMap()), 
00261         std::runtime_error, "Tpetra::Vector::elementwiseMultiply(x,y): Vectors must have compatible Maps.");
00262     TEST_FOR_EXCEPTION( !getMap().isCompatible(y.getMap()), 
00263         std::runtime_error, "Tpetra::Vector::elementwiseMultiply(x,y): Vectors must have compatible Maps.");
00264     // do work
00265     // update flops counter
00266     TEST_FOR_EXCEPT(true);
00267   }
00268 
00269 
00270   template<typename Ordinal, typename Scalar>
00271   void Vector<Ordinal,Scalar>::elementwiseReciprocalMultiply(
00272       Scalar scalarXY, const Vector<Ordinal, Scalar> &x, const Vector<Ordinal, Scalar> &y, 
00273       const Scalar &scalarThis) 
00274   {
00275     TEST_FOR_EXCEPTION( !getMap().isCompatible(x.getMap()), 
00276         std::runtime_error, "Tpetra::Vector::elementwiseReciprocalMultiply(x,y): Vectors must have compatible Maps.");
00277     TEST_FOR_EXCEPTION( !getMap().isCompatible(y.getMap()), 
00278         std::runtime_error, "Tpetra::Vector::elementwiseReciprocalMultiply(x,y): Vectors must have compatible Maps.");
00279     // do work
00280     // update flops counter
00281     TEST_FOR_EXCEPT(true);
00282   }
00283 
00284 
00285   template<typename Ordinal, typename Scalar>
00286   const Scalar & Vector<Ordinal,Scalar>::getSeed() const {
00287     return(VectorData_->seed_);
00288   }
00289 
00290 
00291   template<typename Ordinal, typename Scalar>
00292   void Vector<Ordinal,Scalar>::setSeed(const Scalar &seed) {
00293     VectorData_->seed_ = seed;
00294   }
00295 
00296 
00297 
00298   template<typename Ordinal, typename Scalar>
00299   Scalar& Vector<Ordinal,Scalar>::operator[](Ordinal index) {
00300     return VectorData_->scalarArray_[index];
00301   }
00302 
00303 
00304   template<typename Ordinal, typename Scalar>
00305   const Scalar & Vector<Ordinal,Scalar>::operator[](Ordinal index) const {
00306     return VectorData_->scalarArray_[index];
00307   }
00308 
00309 
00310   template<typename Ordinal, typename Scalar>
00311   Ordinal Vector<Ordinal,Scalar>::getNumMyEntries() const {
00312     return VectorData_->map_.getNumMyEntries();
00313   }
00314 
00315 
00316   template<typename Ordinal, typename Scalar>
00317   Ordinal Vector<Ordinal,Scalar>::getNumGlobalEntries() const {
00318     return VectorData_->map_.getNumGlobalEntries();
00319   }
00320 
00321 
00322   template<typename Ordinal, typename Scalar>
00323   void Vector<Ordinal,Scalar>::print(std::ostream& os) const {
00324     TEST_FOR_EXCEPT(true);
00325   }
00326 
00327   template<typename Ordinal, typename Scalar>
00328   void Vector<Ordinal,Scalar>::printValues(std::ostream& os) const {
00329     TEST_FOR_EXCEPT(true);
00330   }
00331 
00332 
00333   template<typename Ordinal, typename Scalar>
00334   const Map<Ordinal> & Vector<Ordinal,Scalar>::getMap() const {
00335     return VectorData_->map_;
00336   }
00337 
00338 
00339   template<typename Ordinal, typename Scalar>
00340   Vector<Ordinal,Scalar> & Vector<Ordinal,Scalar>::operator=(const Vector<Ordinal,Scalar> &source) 
00341   {
00342     VectorData_ = source.VectorData_;
00343     return *this;
00344   }
00345 
00346 
00347   template<typename Ordinal, typename Scalar>
00348   Teuchos::ArrayView<Scalar> Vector<Ordinal,Scalar>::scalarPointer() 
00349   {
00350     return VectorData_->values_;
00351   }
00352 
00353   template<typename Ordinal, typename Scalar>
00354   Teuchos::ArrayView<const Scalar> Vector<Ordinal,Scalar>::scalarPointer() const
00355   {
00356     return VectorData_->values_;
00357   }
00358 
00359   template<typename Ordinal, typename Scalar>
00360   bool Vector<Ordinal,Scalar>::checkSizes(const DistObject<Ordinal,Scalar> &sourceObj) {
00361     // first check that sourceObj is actually a Vector, and not some other kind of DistObject
00362     /*try {
00363       Vector<Ordinal, Scalar> const& sourceVector = dynamic_cast<Vector<Ordinal, Scalar> const&>(sourceObj);
00364       }
00365       catch(std::bad_cast bc) {
00366       return(false);
00367       }*/
00368 
00369     // ???
00370     TEST_FOR_EXCEPT(true);
00371     return true;
00372   }
00373 
00374   template<typename Ordinal, typename Scalar>
00375   void Vector<Ordinal,Scalar>::copyAndPermute(
00376       const DistObject<Ordinal,Scalar> &sourceObj,
00377       Ordinal numSameIDs, Ordinal numPermuteIDs,
00378       const Teuchos::ArrayView<const Ordinal> &permuteToLIDs, 
00379       const Teuchos::ArrayView<const Ordinal> &permuteFromLIDs) 
00380   {
00381     // cast sourceObj to a Tpetra::Vector so we can actually do something with it
00382     // const Vector<Ordinal,Scalar> &sourceVector = dynamic_cast<const Vector<Ordinal,Scalar> &>(sourceObj);
00383 
00384     /*
00385     const std::vector<Scalar> &sourceArray = sourceVector.scalarArray();
00386     std::vector<Scalar> &destArray = scalarArray();
00387 
00388     // the first numImportIDs GIDs are the same between source and target,
00389     // we can just copy them
00390     for(Ordinal i = Teuchos::OrdinalTraits<Ordinal>::zero(); i < numSameIDs; i++)
00391       destArray[i] = sourceArray[i];
00392 
00393     // next, do permutations
00394     for(Ordinal i = Teuchos::OrdinalTraits<Ordinal>::zero(); i < numPermuteIDs; i++)
00395       destArray[permuteToLIDs[i]] = sourceArray[permuteFromLIDs[i]];
00396 
00397      */
00398     TEST_FOR_EXCEPT(true);
00399   }
00400 
00401   template<typename Ordinal, typename Scalar>
00402   void Vector<Ordinal,Scalar>::packAndPrepare(
00403       const DistObject<Ordinal,Scalar> &sourceObj,
00404       Ordinal numExportIDs,
00405       const Teuchos::ArrayView<const Ordinal> &exportLIDs, 
00406       const Teuchos::ArrayView<Scalar> &exports,
00407       Ordinal &packetSize,
00408       Distributor<Ordinal> &distor) 
00409   {
00410     /*
00411     // cast sourceObj to a Tpetra::Vector so we can actually do something with it
00412     Vector<Ordinal, Scalar> const& sourceVector = dynamic_cast<Vector<Ordinal, Scalar> const&>(sourceObj);
00413 
00414     // For a vector, we only need to send the value
00415     exports.clear();
00416     for(Ordinal i = Teuchos::OrdinalTraits<Ordinal>::zero(); i < numExportIDs; i++)
00417       exports.push_back(sourceVector[exportLIDs[i]]);
00418 
00419     // packetSize = 1
00420     packetSize = Teuchos::OrdinalTraits<Ordinal>::one();
00421     */
00422     TEST_FOR_EXCEPT(true);
00423   }
00424 
00425   template<typename Ordinal, typename Scalar>
00426   void Vector<Ordinal,Scalar>::unpackAndCombine(
00427       Ordinal numImportIDs,
00428       const Teuchos::ArrayView<const Ordinal> &importLIDs,
00429       const Teuchos::ArrayView<const Scalar> &imports,
00430       Distributor<Ordinal> &distor,
00431       CombineMode CM) 
00432   {
00433     /*
00434     if(CM == Insert || CM == Replace) {
00435       // copy values from scalarExports
00436       for(Ordinal i = Teuchos::OrdinalTraits<Ordinal>::zero(); i < numImportIDs; i++)
00437         scalarArray().at(importLIDs[i]) = imports[i];
00438     }
00439     else if(CM == Add) {
00440       // sum values from scalarExports
00441       submitEntries(numImportIDs, &importLIDs.front(), &imports.front());
00442     }
00443     else
00444       throw Teuchos::Object::reportError("Unknown CombineMode", -99);
00445       */
00446 
00447     TEST_FOR_EXCEPT(true);
00448   }
00449 
00450 } // namespace Tpetra
00451 
00452 #endif // TPETRA_VECTOR_HPP
Tpetra_Vector.hpp
