00001
00002
00003
00004
00005
00006
00007
00008
00009
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
00043 template <typename number> class Vector;
00044
00045
00050 namespace TrilinosWrappers
00051 {
00052
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
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
01060
01061
01062
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
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
01214
01215
01216
01217
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
01232
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
01291
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
01325
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
01341
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
01358
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
01515
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
01531
01532
01533
01534
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
01545
01546
01547
01548
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
01749
01750
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
01781 if (vector->Map().SameAs(v.vector->Map())==false)
01782 {
01783 *vector = *v.vector;
01784 *this *= a;
01785 }
01786 else
01787 {
01788
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
01813 if (vector->Map().SameAs(v.vector->Map())==false)
01814 {
01815 *vector = *v.vector;
01816 sadd(a, b, w);
01817 }
01818 else
01819 {
01820
01821
01822
01823
01824
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
01891
01892 #endif
01893
01894