include/deal.II/lac/trilinos_vector_base.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: trilinos_vector_base.h 25526 2012-05-21 12:08:58Z bangerth @f$
00003 //
00004 //    Copyright (C) 2008, 2009, 2010, 2011, 2012 by the deal.II authors
00005 //
00006 //    This file is subject to QPL and may not be  distributed
00007 //    without copyright and license information. Please refer
00008 //    to the file deal.II/doc/license.html for the  text  and
00009 //    further information on this license.
00010 //
00011 //---------------------------------------------------------------------------
00012 #ifndef __deal2__trilinos_vector_base_h
00013 #define __deal2__trilinos_vector_base_h
00014 
00015 
00016 #include <deal.II/base/config.h>
00017 
00018 #ifdef DEAL_II_USE_TRILINOS
00019 
00020 #include <deal.II/base/utilities.h>
00021 #  include <deal.II/base/std_cxx1x/shared_ptr.h>
00022 #  include <deal.II/base/subscriptor.h>
00023 #  include <deal.II/lac/exceptions.h>
00024 #  include <deal.II/lac/vector.h>
00025 
00026 #  include <vector>
00027 #  include <utility>
00028 #  include <memory>
00029 
00030 #  define TrilinosScalar double
00031 #  include "Epetra_ConfigDefs.h"
00032 #  ifdef DEAL_II_COMPILER_SUPPORTS_MPI // only if MPI is installed
00033 #    include "mpi.h"
00034 #    include "Epetra_MpiComm.h"
00035 #  else
00036 #    include "Epetra_SerialComm.h"
00037 #  endif
00038 #  include "Epetra_FEVector.h"
00039 
00040 DEAL_II_NAMESPACE_OPEN
00041 
00042                                 // forward declaration
00043 template <typename number> class Vector;
00044 
00045 
00050 namespace TrilinosWrappers
00051 {
00052                                 // forward declaration
00053   class VectorBase;
00054 
00055 
00066   namespace internal
00067   {
00084     class VectorReference
00085     {
00086       private:
00093         VectorReference (VectorBase        &vector,
00094                          const unsigned int index);
00095 
00096       public:
00121         const VectorReference &
00122           operator = (const VectorReference &r) const;
00123 
00128         const VectorReference &
00129           operator = (const VectorReference &r);
00130 
00135         const VectorReference &
00136           operator = (const TrilinosScalar &s) const;
00137 
00143         const VectorReference &
00144           operator += (const TrilinosScalar &s) const;
00145 
00151         const VectorReference &
00152           operator -= (const TrilinosScalar &s) const;
00153 
00159         const VectorReference &
00160           operator *= (const TrilinosScalar &s) const;
00161 
00167         const VectorReference &
00168           operator /= (const TrilinosScalar &s) const;
00169 
00176         operator TrilinosScalar () const;
00177 
00181         DeclException1 (ExcTrilinosError,
00182                         int,
00183                         << "An error with error number " << arg1
00184                         << " occured while calling a Trilinos function");
00185 
00189         DeclException3 (ExcAccessToNonLocalElement,
00190                         int, int, int,
00191                         << "You tried to access element " << arg1
00192                         << " of a distributed vector, but only elements "
00193                         << arg2 << " through " << arg3
00194                         << " are stored locally and can be accessed.");
00195 
00196       private:
00201         VectorBase   &vector;
00202 
00207         const unsigned int  index;
00208 
00215         friend class ::TrilinosWrappers::VectorBase;
00216     };
00217   }
00255   class VectorBase : public Subscriptor
00256   {
00257     public:
00266       typedef TrilinosScalar            value_type;
00267       typedef TrilinosScalar            real_type;
00268       typedef std::size_t               size_type;
00269       typedef internal::VectorReference reference;
00270       typedef const internal::VectorReference const_reference;
00271 
00276 
00287       VectorBase ();
00288 
00295       VectorBase (const VectorBase &v);
00296 
00300       virtual ~VectorBase ();
00301 
00308       void clear ();
00309 
00317       void reinit (const VectorBase &v,
00318                    const bool        fast = false);
00319 
00349       void compress (const Epetra_CombineMode last_action = Zero);
00350 
00358       bool is_compressed () const;
00359 
00383       VectorBase &
00384         operator = (const TrilinosScalar s);
00385 
00394       VectorBase &
00395         operator = (const VectorBase &v);
00396 
00414       template <typename Number>
00415       VectorBase &
00416         operator = (const ::Vector<Number> &v);
00417 
00427       bool operator == (const VectorBase &v) const;
00428 
00438       bool operator != (const VectorBase &v) const;
00439 
00444       unsigned int size () const;
00445 
00460       unsigned int local_size () const;
00461 
00477       std::pair<unsigned int, unsigned int> local_range () const;
00478 
00484       bool in_local_range (const unsigned int index) const;
00485 
00490       bool has_ghost_elements() const;
00491 
00498       TrilinosScalar operator * (const VectorBase &vec) const;
00499 
00504       real_type norm_sqr () const;
00505 
00510       TrilinosScalar mean_value () const;
00511 
00516       TrilinosScalar minimal_value () const;
00517 
00522       real_type l1_norm () const;
00523 
00529       real_type l2_norm () const;
00530 
00538       real_type lp_norm (const TrilinosScalar p) const;
00539 
00544       real_type linfty_norm () const;
00545 
00556       bool all_zero () const;
00557 
00568       bool is_non_negative () const;
00570 
00571 
00576 
00581       reference
00582         operator () (const unsigned int index);
00583 
00589       TrilinosScalar
00590         operator () (const unsigned int index) const;
00591 
00598       reference
00599         operator [] (const unsigned int index);
00600 
00608       TrilinosScalar
00609         operator [] (const unsigned int index) const;
00610 
00622       TrilinosScalar el (const unsigned int index) const;
00623 
00635       void set (const std::vector<unsigned int>    &indices,
00636                 const std::vector<TrilinosScalar>  &values);
00637 
00645       void set (const std::vector<unsigned int>        &indices,
00646                 const ::Vector<TrilinosScalar> &values);
00648 
00649 
00654 
00665       void set (const unsigned int    n_elements,
00666                 const unsigned int   *indices,
00667                 const TrilinosScalar *values);
00668 
00677       void add (const std::vector<unsigned int>   &indices,
00678                 const std::vector<TrilinosScalar> &values);
00679 
00687       void add (const std::vector<unsigned int>        &indices,
00688                 const ::Vector<TrilinosScalar> &values);
00689 
00699       void add (const unsigned int    n_elements,
00700                 const unsigned int   *indices,
00701                 const TrilinosScalar *values);
00702 
00707       VectorBase & operator *= (const TrilinosScalar factor);
00708 
00713       VectorBase & operator /= (const TrilinosScalar factor);
00714 
00719       VectorBase & operator += (const VectorBase &V);
00720 
00725       VectorBase & operator -= (const VectorBase &V);
00726 
00732       void add (const TrilinosScalar s);
00733 
00744       void add (const VectorBase &V,
00745                 const bool        allow_different_maps = false);
00746 
00752       void add (const TrilinosScalar  a,
00753                 const VectorBase     &V);
00754 
00760       void add (const TrilinosScalar  a,
00761                 const VectorBase     &V,
00762                 const TrilinosScalar  b,
00763                 const VectorBase     &W);
00764 
00770       void sadd (const TrilinosScalar  s,
00771                  const VectorBase     &V);
00772 
00778       void sadd (const TrilinosScalar  s,
00779                  const TrilinosScalar  a,
00780                  const VectorBase     &V);
00781 
00785       void sadd (const TrilinosScalar  s,
00786                  const TrilinosScalar  a,
00787                  const VectorBase     &V,
00788                  const TrilinosScalar  b,
00789                  const VectorBase     &W);
00790 
00796       void sadd (const TrilinosScalar  s,
00797                  const TrilinosScalar  a,
00798                  const VectorBase     &V,
00799                  const TrilinosScalar  b,
00800                  const VectorBase     &W,
00801                  const TrilinosScalar  c,
00802                  const VectorBase     &X);
00803 
00813       void scale (const VectorBase &scaling_factors);
00814 
00819       void equ (const TrilinosScalar  a,
00820                 const VectorBase     &V);
00821 
00826       void equ (const TrilinosScalar  a,
00827                 const VectorBase     &V,
00828                 const TrilinosScalar  b,
00829                 const VectorBase     &W);
00830 
00848       void ratio (const VectorBase &a,
00849                   const VectorBase &b);
00851 
00852 
00857 
00863       const Epetra_MultiVector & trilinos_vector () const;
00864 
00870       Epetra_FEVector & trilinos_vector ();
00871 
00878       const Epetra_Map & vector_partitioner () const;
00879 
00886       void print (const char* format = 0) const;
00887 
00901       void print (std::ostream       &out,
00902                   const unsigned int  precision  = 3,
00903                   const bool          scientific = true,
00904                   const bool          across     = true) const;
00905 
00932       void swap (VectorBase &v);
00933 
00938       std::size_t memory_consumption () const;
00940 
00944       DeclException0 (ExcGhostsPresent);
00945 
00949       DeclException0 (ExcDifferentParallelPartitioning);
00950 
00954       DeclException1 (ExcTrilinosError,
00955                       int,
00956                       << "An error with error number " << arg1
00957                       << " occured while calling a Trilinos function");
00958 
00962       DeclException3 (ExcAccessToNonlocalElement,
00963                       int, int, int,
00964                       << "You tried to access element " << arg1
00965                       << " of a distributed vector, but only entries "
00966                       << arg2 << " through " << arg3
00967                       << " are stored locally and can be accessed.");
00968 
00969 
00970     private:
00996       Epetra_CombineMode last_action;
00997 
01003       bool compressed;
01004 
01010       std_cxx1x::shared_ptr<Epetra_FEVector> vector;
01011 
01012 
01017       friend class internal::VectorReference;
01018       friend class Vector;
01019       friend class MPI::Vector;
01020   };
01021 
01022 
01023 
01024 
01025 // ------------------- inline and template functions --------------
01026 
01035   inline
01036   void swap (VectorBase &u, VectorBase &v)
01037   {
01038     u.swap (v);
01039   }
01040 
01041 
01042 #ifndef DOXYGEN
01043 
01044   namespace internal
01045   {
01046     inline
01047     VectorReference::VectorReference (VectorBase        &vector,
01048                                       const unsigned int index)
01049                     :
01050                     vector (vector),
01051                     index (index)
01052     {}
01053 
01054 
01055     inline
01056     const VectorReference &
01057     VectorReference::operator = (const VectorReference &r) const
01058     {
01059                                         // as explained in the class
01060                                         // documentation, this is not the copy
01061                                         // operator. so simply pass on to the
01062                                         // "correct" assignment operator
01063       *this = static_cast<TrilinosScalar> (r);
01064 
01065       return *this;
01066     }
01067 
01068 
01069 
01070     inline
01071     const VectorReference &
01072     VectorReference::operator = (const VectorReference &r)
01073     {
01074                                         // as above
01075       *this = static_cast<TrilinosScalar> (r);
01076 
01077       return *this;
01078     }
01079 
01080 
01081     inline
01082     const VectorReference &
01083     VectorReference::operator = (const TrilinosScalar &value) const
01084     {
01085       vector.set (1, &index, &value);
01086       return *this;
01087     }
01088 
01089 
01090 
01091     inline
01092     const VectorReference &
01093     VectorReference::operator += (const TrilinosScalar &value) const
01094     {
01095       vector.add (1, &index, &value);
01096       return *this;
01097     }
01098 
01099 
01100 
01101     inline
01102     const VectorReference &
01103     VectorReference::operator -= (const TrilinosScalar &value) const
01104     {
01105       TrilinosScalar new_value = -value;
01106       vector.add (1, &index, &new_value);
01107       return *this;
01108     }
01109 
01110 
01111 
01112     inline
01113     const VectorReference &
01114     VectorReference::operator *= (const TrilinosScalar &value) const
01115     {
01116       TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
01117       vector.set (1, &index, &new_value);
01118       return *this;
01119     }
01120 
01121 
01122 
01123     inline
01124     const VectorReference &
01125     VectorReference::operator /= (const TrilinosScalar &value) const
01126     {
01127       TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
01128       vector.set (1, &index, &new_value);
01129       return *this;
01130     }
01131   }
01132 
01133 
01134 
01135   inline
01136   bool
01137   VectorBase::is_compressed () const
01138   {
01139     return compressed;
01140   }
01141 
01142 
01143 
01144   inline
01145   bool
01146   VectorBase::in_local_range (const unsigned int index) const
01147   {
01148     std::pair<unsigned int, unsigned int> range = local_range();
01149 
01150     return ((index >= range.first) && (index <  range.second));
01151   }
01152 
01153 
01154 
01155   inline
01156   bool
01157   VectorBase::has_ghost_elements() const
01158   {
01159     return vector->Map().UniqueGIDs()==false;
01160   }
01161 
01162 
01163 
01164   inline
01165   internal::VectorReference
01166   VectorBase::operator () (const unsigned int index)
01167   {
01168     return internal::VectorReference (*this, index);
01169   }
01170 
01171 
01172 
01173   inline
01174   internal::VectorReference
01175   VectorBase::operator [] (const unsigned int index)
01176   {
01177     return operator() (index);
01178   }
01179 
01180 
01181   inline
01182   TrilinosScalar
01183   VectorBase::operator [] (const unsigned int index) const
01184   {
01185     return operator() (index);
01186   }
01187 
01188 
01189 
01190   inline
01191   void
01192   VectorBase::reinit (const VectorBase &v,
01193                       const bool        fast)
01194   {
01195     Assert (vector.get() != 0,
01196             ExcMessage("Vector has not been constructed properly."));
01197 
01198     if (fast == false ||
01199         vector_partitioner().SameAs(v.vector_partitioner())==false)
01200       vector.reset (new Epetra_FEVector(*v.vector));
01201   }
01202 
01203 
01204 
01205   inline
01206   void
01207   VectorBase::compress (const Epetra_CombineMode given_last_action)
01208   {
01209     Epetra_CombineMode mode = (last_action != Zero) ?
01210                               last_action : given_last_action;
01211 #ifdef DEBUG
01212 #  ifdef DEAL_II_COMPILER_SUPPORTS_MPI
01213                                      // check that every process has decided
01214                                      // to use the same mode. This will
01215                                      // otherwise result in undefined
01216                                      // behaviour in the call to
01217                                      // GlobalAssemble().
01218     double double_mode = mode;
01219     Utilities::MPI::MinMaxAvg result
01220       = Utilities::MPI::min_max_avg (double_mode,
01221                                      dynamic_cast<const Epetra_MpiComm*>
01222                                      (&vector_partitioner().Comm())->GetMpiComm());
01223     Assert(result.max-result.min<1e-5,
01224            ExcMessage ("Not all processors agree whether the last operation on "
01225                        "this vector was an addition or a set operation. This will "
01226                        "prevent the compress() operation from succeeding."));
01227 
01228 #  endif
01229 #endif
01230 
01231                                  // Now pass over the information about
01232                                  // what we did last to the vector.
01233     const int ierr = vector->GlobalAssemble(mode);
01234     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01235     last_action = Zero;
01236 
01237     compressed = true;
01238   }
01239 
01240 
01241 
01242   inline
01243   VectorBase &
01244   VectorBase::operator = (const TrilinosScalar s)
01245   {
01246 
01247     Assert (numbers::is_finite(s), ExcNumberNotFinite());
01248 
01249     const int ierr = vector->PutScalar(s);
01250 
01251     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01252 
01253     return *this;
01254   }
01255 
01256 
01257 
01258   inline
01259   void
01260   VectorBase::set (const std::vector<unsigned int>    &indices,
01261                    const std::vector<TrilinosScalar>  &values)
01262   {
01263     Assert (indices.size() == values.size(),
01264             ExcDimensionMismatch(indices.size(),values.size()));
01265 
01266     set (indices.size(), &indices[0], &values[0]);
01267   }
01268 
01269 
01270 
01271   inline
01272   void
01273   VectorBase::set (const std::vector<unsigned int>        &indices,
01274                    const ::Vector<TrilinosScalar> &values)
01275   {
01276     Assert (indices.size() == values.size(),
01277             ExcDimensionMismatch(indices.size(),values.size()));
01278 
01279     set (indices.size(), &indices[0], values.begin());
01280   }
01281 
01282 
01283 
01284   inline
01285   void
01286   VectorBase::set (const unsigned int    n_elements,
01287                    const unsigned int   *indices,
01288                    const TrilinosScalar *values)
01289   {
01290                                      // if we have ghost values, do not allow
01291                                      // writing to this vector at all.
01292     Assert (!has_ghost_elements(), ExcGhostsPresent());
01293 
01294     if (last_action == Add)
01295       vector->GlobalAssemble(Add);
01296 
01297     if (last_action != Insert)
01298       last_action = Insert;
01299 
01300     for (unsigned int i=0; i<n_elements; ++i)
01301       {
01302         const unsigned int row = indices[i];
01303         const int local_row = vector->Map().LID(indices[i]);
01304         if (local_row == -1)
01305           {
01306             const int ierr = vector->ReplaceGlobalValues (1,
01307                                                           (const int*)(&row),
01308                                                           &values[i]);
01309             AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01310             compressed = false;
01311           }
01312         else
01313           (*vector)[0][local_row] = values[i];
01314       }
01315   }
01316 
01317 
01318 
01319   inline
01320   void
01321   VectorBase::add (const std::vector<unsigned int>    &indices,
01322                    const std::vector<TrilinosScalar>  &values)
01323   {
01324                                      // if we have ghost values, do not allow
01325                                      // writing to this vector at all.
01326     Assert (!has_ghost_elements(), ExcGhostsPresent());
01327     Assert (indices.size() == values.size(),
01328             ExcDimensionMismatch(indices.size(),values.size()));
01329 
01330     add (indices.size(), &indices[0], &values[0]);
01331   }
01332 
01333 
01334 
01335   inline
01336   void
01337   VectorBase::add (const std::vector<unsigned int>        &indices,
01338                    const ::Vector<TrilinosScalar> &values)
01339   {
01340                                      // if we have ghost values, do not allow
01341                                      // writing to this vector at all.
01342     Assert (!has_ghost_elements(), ExcGhostsPresent());
01343     Assert (indices.size() == values.size(),
01344             ExcDimensionMismatch(indices.size(),values.size()));
01345 
01346     add (indices.size(), &indices[0], values.begin());
01347   }
01348 
01349 
01350 
01351   inline
01352   void
01353   VectorBase::add (const unsigned int    n_elements,
01354                    const unsigned int   *indices,
01355                    const TrilinosScalar *values)
01356   {
01357                                      // if we have ghost values, do not allow
01358                                      // writing to this vector at all.
01359     Assert (!has_ghost_elements(), ExcGhostsPresent());
01360 
01361     if (last_action != Add)
01362       {
01363         if (last_action == Insert)
01364           vector->GlobalAssemble(Insert);
01365         last_action = Add;
01366       }
01367 
01368     for (unsigned int i=0; i<n_elements; ++i)
01369       {
01370         const unsigned int row = indices[i];
01371         const int local_row = vector->Map().LID(row);
01372         if (local_row == -1)
01373           {
01374             const int ierr = vector->SumIntoGlobalValues (1,
01375                                                           (const int*)(&row),
01376                                                           &values[i]);
01377             AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01378             compressed = false;
01379           }
01380         else
01381           (*vector)[0][local_row] += values[i];
01382       }
01383   }
01384 
01385 
01386 
01387   inline
01388   unsigned int
01389   VectorBase::size () const
01390   {
01391     return (unsigned int) (vector->Map().MaxAllGID() + 1 -
01392                            vector->Map().MinAllGID());
01393   }
01394 
01395 
01396 
01397   inline
01398   unsigned int
01399   VectorBase::local_size () const
01400   {
01401     return (unsigned int) vector->Map().NumMyElements();
01402   }
01403 
01404 
01405 
01406   inline
01407   std::pair<unsigned int, unsigned int>
01408   VectorBase::local_range () const
01409   {
01410     int begin, end;
01411     begin = vector->Map().MinMyGID();
01412     end = vector->Map().MaxMyGID()+1;
01413     return std::make_pair (begin, end);
01414   }
01415 
01416 
01417 
01418   inline
01419   TrilinosScalar
01420   VectorBase::operator * (const VectorBase &vec) const
01421   {
01422     Assert (vector->Map().SameAs(vec.vector->Map()),
01423             ExcDifferentParallelPartitioning());
01424     Assert (!has_ghost_elements(), ExcGhostsPresent());
01425 
01426     TrilinosScalar result;
01427 
01428     const int ierr = vector->Dot(*(vec.vector), &result);
01429     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01430 
01431     return result;
01432   }
01433 
01434 
01435 
01436   inline
01437   VectorBase::real_type
01438   VectorBase::norm_sqr () const
01439   {
01440     const TrilinosScalar d = l2_norm();
01441     return d*d;
01442   }
01443 
01444 
01445 
01446   inline
01447   TrilinosScalar
01448   VectorBase::mean_value () const
01449   {
01450     Assert (!has_ghost_elements(), ExcGhostsPresent());
01451 
01452     TrilinosScalar mean;
01453     const int ierr = vector->MeanValue (&mean);
01454     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01455 
01456     return mean;
01457   }
01458 
01459 
01460 
01461   inline
01462   TrilinosScalar
01463   VectorBase::minimal_value () const
01464   {
01465     TrilinosScalar min_value;
01466     const int ierr = vector->MinValue (&min_value);
01467     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01468 
01469     return min_value;
01470   }
01471 
01472 
01473 
01474   inline
01475   VectorBase::real_type
01476   VectorBase::l1_norm () const
01477   {
01478     Assert (!has_ghost_elements(), ExcGhostsPresent());
01479 
01480     TrilinosScalar d;
01481     const int ierr = vector->Norm1 (&d);
01482     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01483 
01484     return d;
01485   }
01486 
01487 
01488 
01489   inline
01490   VectorBase::real_type
01491   VectorBase::l2_norm () const
01492   {
01493     Assert (!has_ghost_elements(), ExcGhostsPresent());
01494 
01495     TrilinosScalar d;
01496     const int ierr = vector->Norm2 (&d);
01497     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01498 
01499     return d;
01500   }
01501 
01502 
01503 
01504   inline
01505   VectorBase::real_type
01506   VectorBase::lp_norm (const TrilinosScalar p) const
01507   {
01508     Assert (!has_ghost_elements(), ExcGhostsPresent());
01509 
01510     TrilinosScalar norm = 0;
01511     TrilinosScalar sum=0;
01512     const unsigned int n_local = local_size();
01513 
01514                                         // loop over all the elements because
01515                                         // Trilinos does not support lp norms
01516     for (unsigned int i=0; i<n_local; i++)
01517       sum += std::pow(std::fabs((*vector)[0][i]), p);
01518 
01519     norm = std::pow(sum, static_cast<TrilinosScalar>(1./p));
01520 
01521     return norm;
01522   }
01523 
01524 
01525 
01526   inline
01527   VectorBase::real_type
01528   VectorBase::linfty_norm () const
01529   {
01530                                      // while we disallow the other
01531                                      // norm operations on ghosted
01532                                      // vectors, this particular norm
01533                                      // is safe to run even in the
01534                                      // presence of ghost elements
01535     TrilinosScalar d;
01536     const int ierr = vector->NormInf (&d);
01537     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01538 
01539     return d;
01540   }
01541 
01542 
01543 
01544                                    // inline also scalar products, vector
01545                                    // additions etc. since they are all
01546                                    // representable by a single Trilinos
01547                                    // call. This reduces the overhead of the
01548                                    // wrapper class.
01549   inline
01550   VectorBase &
01551   VectorBase::operator *= (const TrilinosScalar a)
01552   {
01553     Assert (numbers::is_finite(a), ExcNumberNotFinite());
01554 
01555     const int ierr = vector->Scale(a);
01556     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01557 
01558     return *this;
01559   }
01560 
01561 
01562 
01563   inline
01564   VectorBase &
01565   VectorBase::operator /= (const TrilinosScalar a)
01566   {
01567     Assert (numbers::is_finite(a), ExcNumberNotFinite());
01568 
01569     const TrilinosScalar factor = 1./a;
01570 
01571     Assert (numbers::is_finite(factor), ExcNumberNotFinite());
01572 
01573     const int ierr = vector->Scale(factor);
01574     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01575 
01576     return *this;
01577   }
01578 
01579 
01580 
01581   inline
01582   VectorBase &
01583   VectorBase::operator += (const VectorBase &v)
01584   {
01585     Assert (size() == v.size(),
01586             ExcDimensionMismatch(size(), v.size()));
01587     Assert (vector->Map().SameAs(v.vector->Map()),
01588             ExcDifferentParallelPartitioning());
01589 
01590     const int ierr = vector->Update (1.0, *(v.vector), 1.0);
01591     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01592 
01593     return *this;
01594   }
01595 
01596 
01597 
01598   inline
01599   VectorBase &
01600   VectorBase::operator -= (const VectorBase &v)
01601   {
01602     Assert (size() == v.size(),
01603             ExcDimensionMismatch(size(), v.size()));
01604     Assert (vector->Map().SameAs(v.vector->Map()),
01605             ExcDifferentParallelPartitioning());
01606 
01607     const int ierr = vector->Update (-1.0, *(v.vector), 1.0);
01608     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01609 
01610     return *this;
01611   }
01612 
01613 
01614 
01615   inline
01616   void
01617   VectorBase::add (const TrilinosScalar s)
01618   {
01619     Assert (numbers::is_finite(s), ExcNumberNotFinite());
01620 
01621     unsigned int n_local = local_size();
01622     for (unsigned int i=0; i<n_local; i++)
01623       (*vector)[0][i] += s;
01624   }
01625 
01626 
01627 
01628   inline
01629   void
01630   VectorBase::add (const TrilinosScalar  a,
01631                    const VectorBase     &v)
01632   {
01633     Assert (local_size() == v.local_size(),
01634             ExcDimensionMismatch(local_size(), v.local_size()));
01635 
01636     Assert (numbers::is_finite(a), ExcNumberNotFinite());
01637 
01638     const int ierr = vector->Update(a, *(v.vector), 1.);
01639     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01640   }
01641 
01642 
01643 
01644   inline
01645   void
01646   VectorBase::add (const TrilinosScalar  a,
01647                    const VectorBase     &v,
01648                    const TrilinosScalar  b,
01649                    const VectorBase     &w)
01650   {
01651     Assert (local_size() == v.local_size(),
01652             ExcDimensionMismatch(local_size(), v.local_size()));
01653     Assert (local_size() == w.local_size(),
01654             ExcDimensionMismatch(local_size(), w.local_size()));
01655 
01656     Assert (numbers::is_finite(a), ExcNumberNotFinite());
01657     Assert (numbers::is_finite(b), ExcNumberNotFinite());
01658 
01659     const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
01660 
01661     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01662   }
01663 
01664 
01665 
01666   inline
01667   void
01668   VectorBase::sadd (const TrilinosScalar  s,
01669                     const VectorBase     &v)
01670   {
01671     Assert (local_size() == v.local_size(),
01672             ExcDimensionMismatch(local_size(), v.local_size()));
01673 
01674     Assert (numbers::is_finite(s), ExcNumberNotFinite());
01675 
01676     const int ierr = vector->Update(1., *(v.vector), s);
01677 
01678     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01679   }
01680 
01681 
01682 
01683   inline
01684   void
01685   VectorBase::sadd (const TrilinosScalar  s,
01686                     const TrilinosScalar  a,
01687                     const VectorBase     &v)
01688   {
01689     Assert (local_size() == v.local_size(),
01690             ExcDimensionMismatch(local_size(), v.local_size()));
01691 
01692     Assert (numbers::is_finite(s), ExcNumberNotFinite());
01693     Assert (numbers::is_finite(a), ExcNumberNotFinite());
01694 
01695     const int ierr = vector->Update(a, *(v.vector), s);
01696 
01697     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01698   }
01699 
01700 
01701 
01702   inline
01703   void
01704   VectorBase::sadd (const TrilinosScalar  s,
01705                     const TrilinosScalar  a,
01706                     const VectorBase     &v,
01707                     const TrilinosScalar  b,
01708                     const VectorBase     &w)
01709   {
01710     Assert (local_size() == v.local_size(),
01711             ExcDimensionMismatch(local_size(), v.local_size()));
01712     Assert (local_size() == w.local_size(),
01713             ExcDimensionMismatch(local_size(), w.local_size()));
01714 
01715     Assert (numbers::is_finite(s), ExcNumberNotFinite());
01716     Assert (numbers::is_finite(a), ExcNumberNotFinite());
01717     Assert (numbers::is_finite(b), ExcNumberNotFinite());
01718 
01719     const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s);
01720 
01721     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01722   }
01723 
01724 
01725 
01726   inline
01727   void
01728   VectorBase::sadd (const TrilinosScalar  s,
01729                     const TrilinosScalar  a,
01730                     const VectorBase     &v,
01731                     const TrilinosScalar  b,
01732                     const VectorBase     &w,
01733                     const TrilinosScalar  c,
01734                     const VectorBase     &x)
01735   {
01736     Assert (local_size() == v.local_size(),
01737             ExcDimensionMismatch(local_size(), v.local_size()));
01738     Assert (local_size() == w.local_size(),
01739             ExcDimensionMismatch(local_size(), w.local_size()));
01740     Assert (local_size() == x.local_size(),
01741             ExcDimensionMismatch(local_size(), x.local_size()));
01742 
01743     Assert (numbers::is_finite(s), ExcNumberNotFinite());
01744     Assert (numbers::is_finite(a), ExcNumberNotFinite());
01745     Assert (numbers::is_finite(b), ExcNumberNotFinite());
01746     Assert (numbers::is_finite(c), ExcNumberNotFinite());
01747 
01748                                         // Update member can only
01749                                         // input two other vectors so
01750                                         // do it in two steps
01751     const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s);
01752     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01753 
01754     const int jerr = vector->Update(c, *(x.vector), 1.);
01755     Assert (jerr == 0, ExcTrilinosError(jerr));
01756   }
01757 
01758 
01759 
01760   inline
01761   void
01762   VectorBase::scale (const VectorBase &factors)
01763   {
01764     Assert (local_size() == factors.local_size(),
01765             ExcDimensionMismatch(local_size(), factors.local_size()));
01766 
01767     const int ierr = vector->Multiply (1.0, *(factors.vector), *vector, 0.0);
01768     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01769   }
01770 
01771 
01772 
01773   inline
01774   void
01775   VectorBase::equ (const TrilinosScalar  a,
01776                    const VectorBase     &v)
01777   {
01778     Assert (numbers::is_finite(a), ExcNumberNotFinite());
01779 
01780                                    // If we don't have the same map, copy.
01781     if (vector->Map().SameAs(v.vector->Map())==false)
01782       {
01783         *vector = *v.vector;
01784         *this *= a;
01785       }
01786     else
01787       {
01788                                    // Otherwise, just update
01789         int ierr = vector->Update(a, *v.vector, 0.0);
01790         AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01791 
01792         last_action = Zero;
01793       }
01794 
01795   }
01796 
01797 
01798 
01799   inline
01800   void
01801   VectorBase::equ (const TrilinosScalar  a,
01802                    const VectorBase     &v,
01803                    const TrilinosScalar  b,
01804                    const VectorBase     &w)
01805   {
01806     Assert (v.local_size() == w.local_size(),
01807             ExcDimensionMismatch (v.local_size(), w.local_size()));
01808 
01809     Assert (numbers::is_finite(a), ExcNumberNotFinite());
01810     Assert (numbers::is_finite(b), ExcNumberNotFinite());
01811 
01812                                    // If we don't have the same map, copy.
01813      if (vector->Map().SameAs(v.vector->Map())==false)
01814       {
01815         *vector = *v.vector;
01816         sadd(a, b, w);
01817       }
01818     else
01819       {
01820                                    // Otherwise, just update. verify
01821                                    // that *this does not only have
01822                                    // the same map as v (the
01823                                    // if-condition above) but also as
01824                                    // w
01825         Assert (vector->Map().SameAs(w.vector->Map()),
01826                 ExcDifferentParallelPartitioning());
01827         int ierr = vector->Update(a, *v.vector, b, *w.vector, 0.0);
01828         AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01829 
01830         last_action = Zero;
01831       }
01832   }
01833 
01834 
01835 
01836   inline
01837   void
01838   VectorBase::ratio (const VectorBase &v,
01839                      const VectorBase &w)
01840   {
01841     Assert (v.local_size() == w.local_size(),
01842             ExcDimensionMismatch (v.local_size(), w.local_size()));
01843     Assert (local_size() == w.local_size(),
01844             ExcDimensionMismatch (local_size(), w.local_size()));
01845 
01846     const int ierr = vector->ReciprocalMultiply(1.0, *(w.vector),
01847                                                 *(v.vector), 0.0);
01848 
01849     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01850   }
01851 
01852 
01853 
01854   inline
01855   const Epetra_MultiVector &
01856   VectorBase::trilinos_vector () const
01857   {
01858     return static_cast<const Epetra_MultiVector&>(*vector);
01859   }
01860 
01861 
01862 
01863   inline
01864   Epetra_FEVector &
01865   VectorBase::trilinos_vector ()
01866   {
01867     return *vector;
01868   }
01869 
01870 
01871 
01872   inline
01873   const Epetra_Map &
01874   VectorBase::vector_partitioner () const
01875   {
01876     return static_cast<const Epetra_Map&>(vector->Map());
01877   }
01878 
01879 
01880 #endif // DOXYGEN
01881 
01882 }
01883 
01886 DEAL_II_NAMESPACE_CLOSE
01887 
01888 #endif // DEAL_II_USE_TRILINOS
01889 
01890 /*----------------------------   trilinos_vector_base.h     ---------------------------*/
01891 
01892 #endif
01893 /*----------------------------   trilinos_vector_base.h     ---------------------------*/
01894 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Wed May 23 2012 12:03:31 by doxygen 1.7.3