00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__petsc_matrix_base_h
00013 #define __deal2__petsc_matrix_base_h
00014
00015
00016 #include <deal.II/base/config.h>
00017
00018 #ifdef DEAL_II_USE_PETSC
00019
00020 # include <deal.II/base/subscriptor.h>
00021 # include <deal.II/lac/full_matrix.h>
00022 # include <deal.II/lac/exceptions.h>
00023
00024 # include <petscmat.h>
00025 # include <deal.II/base/std_cxx1x/shared_ptr.h>
00026
00027 # include <vector>
00028 # include <cmath>
00029
00030 DEAL_II_NAMESPACE_OPEN
00031
00032 template <typename Matrix> class BlockMatrixBase;
00033
00034
00035 namespace PETScWrappers
00036 {
00037
00038 class VectorBase;
00039 class MatrixBase;
00040
00041 namespace MatrixIterators
00042 {
00057 class const_iterator
00058 {
00059 private:
00063 class Accessor
00064 {
00065 public:
00072 Accessor (const MatrixBase *matrix,
00073 const unsigned int row,
00074 const unsigned int index);
00075
00081 unsigned int row() const;
00082
00088 unsigned int index() const;
00089
00095 unsigned int column() const;
00096
00100 PetscScalar value() const;
00101
00105 DeclException0 (ExcBeyondEndOfMatrix);
00109 DeclException3 (ExcAccessToNonlocalRow,
00110 int, int, int,
00111 << "You tried to access row " << arg1
00112 << " of a distributed matrix, but only rows "
00113 << arg2 << " through " << arg3
00114 << " are stored locally and can be accessed.");
00115
00116 private:
00120 mutable MatrixBase *matrix;
00121
00125 unsigned int a_row;
00126
00130 unsigned int a_index;
00131
00156 std_cxx1x::shared_ptr<const std::vector<unsigned int> > colnum_cache;
00157
00162 std_cxx1x::shared_ptr<const std::vector<PetscScalar> > value_cache;
00163
00171 void visit_present_row ();
00172
00177 friend class const_iterator;
00178 };
00179
00180 public:
00181
00187 const_iterator (const MatrixBase *matrix,
00188 const unsigned int row,
00189 const unsigned int index);
00190
00194 const_iterator& operator++ ();
00195
00199 const_iterator operator++ (int);
00200
00204 const Accessor& operator* () const;
00205
00209 const Accessor* operator-> () const;
00210
00217 bool operator == (const const_iterator&) const;
00221 bool operator != (const const_iterator&) const;
00222
00232 bool operator < (const const_iterator&) const;
00233
00237 DeclException2 (ExcInvalidIndexWithinRow,
00238 int, int,
00239 << "Attempt to access element " << arg2
00240 << " of row " << arg1
00241 << " which doesn't have that many elements.");
00242
00243 private:
00248 Accessor accessor;
00249 };
00250
00251 }
00252
00253
00286 class MatrixBase : public Subscriptor
00287 {
00288 public:
00293 typedef MatrixIterators::const_iterator const_iterator;
00294
00299 typedef PetscScalar value_type;
00300
00304 MatrixBase ();
00305
00310 virtual ~MatrixBase ();
00311
00327 MatrixBase &
00328 operator = (const value_type d);
00335 void clear ();
00336
00352 void set (const unsigned int i,
00353 const unsigned int j,
00354 const PetscScalar value);
00355
00391 void set (const std::vector<unsigned int> &indices,
00392 const FullMatrix<PetscScalar> &full_matrix,
00393 const bool elide_zero_values = false);
00394
00402 void set (const std::vector<unsigned int> &row_indices,
00403 const std::vector<unsigned int> &col_indices,
00404 const FullMatrix<PetscScalar> &full_matrix,
00405 const bool elide_zero_values = false);
00406
00433 void set (const unsigned int row,
00434 const std::vector<unsigned int> &col_indices,
00435 const std::vector<PetscScalar> &values,
00436 const bool elide_zero_values = false);
00437
00464 void set (const unsigned int row,
00465 const unsigned int n_cols,
00466 const unsigned int *col_indices,
00467 const PetscScalar *values,
00468 const bool elide_zero_values = false);
00469
00485 void add (const unsigned int i,
00486 const unsigned int j,
00487 const PetscScalar value);
00488
00526 void add (const std::vector<unsigned int> &indices,
00527 const FullMatrix<PetscScalar> &full_matrix,
00528 const bool elide_zero_values = true);
00529
00537 void add (const std::vector<unsigned int> &row_indices,
00538 const std::vector<unsigned int> &col_indices,
00539 const FullMatrix<PetscScalar> &full_matrix,
00540 const bool elide_zero_values = true);
00541
00569 void add (const unsigned int row,
00570 const std::vector<unsigned int> &col_indices,
00571 const std::vector<PetscScalar> &values,
00572 const bool elide_zero_values = true);
00573
00601 void add (const unsigned int row,
00602 const unsigned int n_cols,
00603 const unsigned int *col_indices,
00604 const PetscScalar *values,
00605 const bool elide_zero_values = true,
00606 const bool col_indices_are_sorted = false);
00607
00638 void clear_row (const unsigned int row,
00639 const PetscScalar new_diag_value = 0);
00640
00654 void clear_rows (const std::vector<unsigned int> &rows,
00655 const PetscScalar new_diag_value = 0);
00656
00674 void compress ();
00675
00692 PetscScalar operator () (const unsigned int i,
00693 const unsigned int j) const;
00694
00706 PetscScalar el (const unsigned int i,
00707 const unsigned int j) const;
00708
00724 PetscScalar diag_element (const unsigned int i) const;
00725
00730 unsigned int m () const;
00731
00736 unsigned int n () const;
00737
00751 unsigned int local_size () const;
00752
00769 std::pair<unsigned int, unsigned int>
00770 local_range () const;
00771
00777 bool in_local_range (const unsigned int index) const;
00778
00785 virtual const MPI_Comm & get_mpi_communicator () const = 0;
00786
00796 unsigned int n_nonzero_elements () const;
00797
00801 unsigned int row_length (const unsigned int row) const;
00802
00815 PetscReal l1_norm () const;
00816
00830 PetscReal linfty_norm () const;
00831
00838 PetscReal frobenius_norm () const;
00839
00840 #if DEAL_II_PETSC_VERSION_GTE(3,1,0)
00841
00847 PetscReal trace () const;
00848 #endif
00849
00854 MatrixBase & operator *= (const PetscScalar factor);
00855
00860 MatrixBase & operator /= (const PetscScalar factor);
00861
00880 void vmult (VectorBase &dst,
00881 const VectorBase &src) const;
00882
00903 void Tvmult (VectorBase &dst,
00904 const VectorBase &src) const;
00905
00926 void vmult_add (VectorBase &dst,
00927 const VectorBase &src) const;
00928
00952 void Tvmult_add (VectorBase &dst,
00953 const VectorBase &src) const;
00954
00991 PetscScalar matrix_norm_square (const VectorBase &v) const;
00992
01015 PetscScalar matrix_scalar_product (const VectorBase &u,
01016 const VectorBase &v) const;
01017
01043 PetscScalar residual (VectorBase &dst,
01044 const VectorBase &x,
01045 const VectorBase &b) const;
01046
01051 const_iterator begin () const;
01052
01056 const_iterator end () const;
01057
01070 const_iterator begin (const unsigned int r) const;
01071
01083 const_iterator end (const unsigned int r) const;
01084
01096 operator Mat () const;
01097
01102 void transpose ();
01103
01111 #if DEAL_II_PETSC_VERSION_LT(3,2,0)
01112 PetscTruth
01113 #else
01114 PetscBool
01115 #endif
01116 is_symmetric (const double tolerance = 1.e-12);
01117
01118 #if DEAL_II_PETSC_VERSION_GTE(2,3,0)
01119
01128 #if DEAL_II_PETSC_VERSION_LT(3,2,0)
01129 PetscTruth
01130 #else
01131 PetscBool
01132 #endif
01133 is_hermitian (const double tolerance = 1.e-12);
01134 #endif
01135
01136
01137
01138
01139
01140
01141
01142 void write_ascii ();
01143
01148 std::size_t memory_consumption() const;
01149
01153 DeclException1 (ExcPETScError,
01154 int,
01155 << "An error with error number " << arg1
01156 << " occured while calling a PETSc function");
01160 DeclException0 (ExcSourceEqualsDestination);
01161
01165 DeclException2 (ExcWrongMode,
01166 int, int,
01167 << "You tried to do a "
01168 << (arg1 == 1 ?
01169 "'set'" :
01170 (arg1 == 2 ?
01171 "'add'" : "???"))
01172 << " operation but the matrix is currently in "
01173 << (arg2 == 1 ?
01174 "'set'" :
01175 (arg2 == 2 ?
01176 "'add'" : "???"))
01177 << " mode. You first have to call 'compress()'.");
01178
01179 protected:
01185 Mat matrix;
01186
01208 struct LastAction
01209 {
01210 enum Values { none, insert, add };
01211 };
01212
01217 LastAction::Values last_action;
01218
01227 void prepare_action(const LastAction::Values new_action);
01228
01251 void prepare_add();
01259 void prepare_set();
01260
01261
01262
01263 private:
01271 #ifdef PETSC_USE_64BIT_INDICES
01272 std::vector<PetscInt> column_indices;
01273 #else
01274 std::vector<int> column_indices;
01275 #endif
01276
01284 std::vector<PetscScalar> column_values;
01285
01286
01292 template <class> friend class ::BlockMatrixBase;
01293 };
01294
01295
01296
01297 #ifndef DOXYGEN
01298
01299
01300
01301 namespace MatrixIterators
01302 {
01303
01304 inline
01305 const_iterator::Accessor::
01306 Accessor (const MatrixBase *matrix,
01307 const unsigned int row,
01308 const unsigned int index)
01309 :
01310 matrix(const_cast<MatrixBase*>(matrix)),
01311 a_row(row),
01312 a_index(index)
01313 {
01314 visit_present_row ();
01315 }
01316
01317
01318 inline
01319 unsigned int
01320 const_iterator::Accessor::row() const
01321 {
01322 Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
01323 return a_row;
01324 }
01325
01326
01327 inline
01328 unsigned int
01329 const_iterator::Accessor::column() const
01330 {
01331 Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
01332 return (*colnum_cache)[a_index];
01333 }
01334
01335
01336 inline
01337 unsigned int
01338 const_iterator::Accessor::index() const
01339 {
01340 Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
01341 return a_index;
01342 }
01343
01344
01345 inline
01346 PetscScalar
01347 const_iterator::Accessor::value() const
01348 {
01349 Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
01350 return (*value_cache)[a_index];
01351 }
01352
01353
01354 inline
01355 const_iterator::
01356 const_iterator(const MatrixBase *matrix,
01357 const unsigned int row,
01358 const unsigned int index)
01359 :
01360 accessor(matrix, row, index)
01361 {}
01362
01363
01364
01365 inline
01366 const_iterator &
01367 const_iterator::operator++ ()
01368 {
01369 Assert (accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
01370
01371 ++accessor.a_index;
01372
01373
01374
01375
01376 if (accessor.a_index >= accessor.colnum_cache->size())
01377 {
01378 accessor.a_index = 0;
01379 ++accessor.a_row;
01380
01381 while ((accessor.a_row < accessor.matrix->m())
01382 &&
01383 (accessor.matrix->row_length(accessor.a_row) == 0))
01384 ++accessor.a_row;
01385
01386 accessor.visit_present_row();
01387 }
01388 return *this;
01389 }
01390
01391
01392 inline
01393 const_iterator
01394 const_iterator::operator++ (int)
01395 {
01396 const const_iterator old_state = *this;
01397 ++(*this);
01398 return old_state;
01399 }
01400
01401
01402 inline
01403 const const_iterator::Accessor &
01404 const_iterator::operator* () const
01405 {
01406 return accessor;
01407 }
01408
01409
01410 inline
01411 const const_iterator::Accessor *
01412 const_iterator::operator-> () const
01413 {
01414 return &accessor;
01415 }
01416
01417
01418 inline
01419 bool
01420 const_iterator::
01421 operator == (const const_iterator& other) const
01422 {
01423 return (accessor.a_row == other.accessor.a_row &&
01424 accessor.a_index == other.accessor.a_index);
01425 }
01426
01427
01428 inline
01429 bool
01430 const_iterator::
01431 operator != (const const_iterator& other) const
01432 {
01433 return ! (*this == other);
01434 }
01435
01436
01437 inline
01438 bool
01439 const_iterator::
01440 operator < (const const_iterator& other) const
01441 {
01442 return (accessor.row() < other.accessor.row() ||
01443 (accessor.row() == other.accessor.row() &&
01444 accessor.index() < other.accessor.index()));
01445 }
01446
01447 }
01448
01449
01450
01451
01452
01453
01454
01455
01456
01457
01458 inline
01459 void
01460 MatrixBase::set (const unsigned int i,
01461 const unsigned int j,
01462 const PetscScalar value)
01463 {
01464 Assert (numbers::is_finite(value), ExcNumberNotFinite());
01465
01466 set (i, 1, &j, &value, false);
01467 }
01468
01469
01470
01471 inline
01472 void
01473 MatrixBase::set (const std::vector<unsigned int> &indices,
01474 const FullMatrix<PetscScalar> &values,
01475 const bool elide_zero_values)
01476 {
01477 Assert (indices.size() == values.m(),
01478 ExcDimensionMismatch(indices.size(), values.m()));
01479 Assert (values.m() == values.n(), ExcNotQuadratic());
01480
01481 for (unsigned int i=0; i<indices.size(); ++i)
01482 set (indices[i], indices.size(), &indices[0], &values(i,0),
01483 elide_zero_values);
01484 }
01485
01486
01487
01488 inline
01489 void
01490 MatrixBase::set (const std::vector<unsigned int> &row_indices,
01491 const std::vector<unsigned int> &col_indices,
01492 const FullMatrix<PetscScalar> &values,
01493 const bool elide_zero_values)
01494 {
01495 Assert (row_indices.size() == values.m(),
01496 ExcDimensionMismatch(row_indices.size(), values.m()));
01497 Assert (col_indices.size() == values.n(),
01498 ExcDimensionMismatch(col_indices.size(), values.n()));
01499
01500 for (unsigned int i=0; i<row_indices.size(); ++i)
01501 set (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
01502 elide_zero_values);
01503 }
01504
01505
01506
01507 inline
01508 void
01509 MatrixBase::set (const unsigned int row,
01510 const std::vector<unsigned int> &col_indices,
01511 const std::vector<PetscScalar> &values,
01512 const bool elide_zero_values)
01513 {
01514 Assert (col_indices.size() == values.size(),
01515 ExcDimensionMismatch(col_indices.size(), values.size()));
01516
01517 set (row, col_indices.size(), &col_indices[0], &values[0],
01518 elide_zero_values);
01519 }
01520
01521
01522
01523 inline
01524 void
01525 MatrixBase::set (const unsigned int row,
01526 const unsigned int n_cols,
01527 const unsigned int *col_indices,
01528 const PetscScalar *values,
01529 const bool elide_zero_values)
01530 {
01531 prepare_action(LastAction::insert);
01532
01533 #ifdef PETSC_USE_64BIT_INDICES
01534 const PetscInt petsc_i = row;
01535 PetscInt * col_index_ptr;
01536 #else
01537 const int petsc_i = row;
01538 int * col_index_ptr;
01539 #endif
01540 PetscScalar const* col_value_ptr;
01541 int n_columns;
01542
01543
01544
01545 #ifndef PETSC_USE_64BIT_INDICES
01546 if (elide_zero_values == false)
01547 {
01548 col_index_ptr = (int*)col_indices;
01549 col_value_ptr = values;
01550 n_columns = n_cols;
01551 }
01552 else
01553 #endif
01554 {
01555
01556
01557 if (column_indices.size() < n_cols)
01558 {
01559 column_indices.resize(n_cols);
01560 column_values.resize(n_cols);
01561 }
01562
01563 n_columns = 0;
01564 for (unsigned int j=0; j<n_cols; ++j)
01565 {
01566 const PetscScalar value = values[j];
01567 Assert (numbers::is_finite(value), ExcNumberNotFinite());
01568 if (value != PetscScalar())
01569 {
01570 column_indices[n_columns] = col_indices[j];
01571 column_values[n_columns] = value;
01572 n_columns++;
01573 }
01574 }
01575 Assert(n_columns <= (int)n_cols, ExcInternalError());
01576
01577 col_index_ptr = &column_indices[0];
01578 col_value_ptr = &column_values[0];
01579 }
01580
01581 const int ierr
01582 = MatSetValues (matrix, 1, &petsc_i, n_columns, col_index_ptr,
01583 col_value_ptr, INSERT_VALUES);
01584 AssertThrow (ierr == 0, ExcPETScError(ierr));
01585 }
01586
01587
01588
01589 inline
01590 void
01591 MatrixBase::add (const unsigned int i,
01592 const unsigned int j,
01593 const PetscScalar value)
01594 {
01595
01596 Assert (numbers::is_finite(value), ExcNumberNotFinite());
01597
01598 if (value == PetscScalar())
01599 {
01600
01601
01602
01603
01604
01605
01606
01607
01608
01609
01610 prepare_action(LastAction::add);
01611
01612 return;
01613 }
01614 else
01615 add (i, 1, &j, &value, false);
01616 }
01617
01618
01619
01620 inline
01621 void
01622 MatrixBase::add (const std::vector<unsigned int> &indices,
01623 const FullMatrix<PetscScalar> &values,
01624 const bool elide_zero_values)
01625 {
01626 Assert (indices.size() == values.m(),
01627 ExcDimensionMismatch(indices.size(), values.m()));
01628 Assert (values.m() == values.n(), ExcNotQuadratic());
01629
01630 for (unsigned int i=0; i<indices.size(); ++i)
01631 add (indices[i], indices.size(), &indices[0], &values(i,0),
01632 elide_zero_values);
01633 }
01634
01635
01636
01637 inline
01638 void
01639 MatrixBase::add (const std::vector<unsigned int> &row_indices,
01640 const std::vector<unsigned int> &col_indices,
01641 const FullMatrix<PetscScalar> &values,
01642 const bool elide_zero_values)
01643 {
01644 Assert (row_indices.size() == values.m(),
01645 ExcDimensionMismatch(row_indices.size(), values.m()));
01646 Assert (col_indices.size() == values.n(),
01647 ExcDimensionMismatch(col_indices.size(), values.n()));
01648
01649 for (unsigned int i=0; i<row_indices.size(); ++i)
01650 add (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
01651 elide_zero_values);
01652 }
01653
01654
01655
01656 inline
01657 void
01658 MatrixBase::add (const unsigned int row,
01659 const std::vector<unsigned int> &col_indices,
01660 const std::vector<PetscScalar> &values,
01661 const bool elide_zero_values)
01662 {
01663 Assert (col_indices.size() == values.size(),
01664 ExcDimensionMismatch(col_indices.size(), values.size()));
01665
01666 add (row, col_indices.size(), &col_indices[0], &values[0],
01667 elide_zero_values);
01668 }
01669
01670
01671
01672 inline
01673 void
01674 MatrixBase::add (const unsigned int row,
01675 const unsigned int n_cols,
01676 const unsigned int *col_indices,
01677 const PetscScalar *values,
01678 const bool elide_zero_values,
01679 const bool )
01680 {
01681 prepare_action(LastAction::add);
01682
01683 #ifdef PETSC_USE_64BIT_INDICES
01684 const PetscInt petsc_i = row;
01685 PetscInt * col_index_ptr;
01686 #else
01687 const int petsc_i = row;
01688 int * col_index_ptr;
01689 #endif
01690 PetscScalar const* col_value_ptr;
01691 int n_columns;
01692
01693
01694
01695 #ifndef PETSC_USE_64BIT_INDICES
01696 if (elide_zero_values == false)
01697 {
01698 col_index_ptr = (int*)col_indices;
01699 col_value_ptr = values;
01700 n_columns = n_cols;
01701 }
01702 else
01703 #endif
01704 {
01705
01706
01707 if (column_indices.size() < n_cols)
01708 {
01709 column_indices.resize(n_cols);
01710 column_values.resize(n_cols);
01711 }
01712
01713 n_columns = 0;
01714 for (unsigned int j=0; j<n_cols; ++j)
01715 {
01716 const PetscScalar value = values[j];
01717 Assert (numbers::is_finite(value), ExcNumberNotFinite());
01718 if (value != PetscScalar())
01719 {
01720 column_indices[n_columns] = col_indices[j];
01721 column_values[n_columns] = value;
01722 n_columns++;
01723 }
01724 }
01725 Assert(n_columns <= (int)n_cols, ExcInternalError());
01726
01727 col_index_ptr = &column_indices[0];
01728 col_value_ptr = &column_values[0];
01729 }
01730
01731 const int ierr
01732 = MatSetValues (matrix, 1, &petsc_i, n_columns, col_index_ptr,
01733 col_value_ptr, ADD_VALUES);
01734 Assert (ierr == 0, ExcPETScError(ierr));
01735 }
01736
01737
01738
01739
01740
01741
01742 inline
01743 PetscScalar
01744 MatrixBase::operator() (const unsigned int i,
01745 const unsigned int j) const
01746 {
01747 return el(i,j);
01748 }
01749
01750
01751
01752 inline
01753 MatrixBase::const_iterator
01754 MatrixBase::begin() const
01755 {
01756 return const_iterator(this, 0, 0);
01757 }
01758
01759
01760 inline
01761 MatrixBase::const_iterator
01762 MatrixBase::end() const
01763 {
01764 return const_iterator(this, m(), 0);
01765 }
01766
01767
01768 inline
01769 MatrixBase::const_iterator
01770 MatrixBase::begin(const unsigned int r) const
01771 {
01772 Assert (r < m(), ExcIndexRange(r, 0, m()));
01773 if (row_length(r) > 0)
01774 return const_iterator(this, r, 0);
01775 else
01776 return end (r);
01777 }
01778
01779
01780 inline
01781 MatrixBase::const_iterator
01782 MatrixBase::end(const unsigned int r) const
01783 {
01784 Assert (r < m(), ExcIndexRange(r, 0, m()));
01785
01786
01787
01788
01789 for (unsigned int i=r+1; i<m(); ++i)
01790 if (row_length(i) > 0)
01791 return const_iterator(this, i, 0);
01792
01793
01794
01795 return end();
01796 }
01797
01798
01799
01800 inline
01801 bool
01802 MatrixBase::in_local_range (const unsigned int index) const
01803 {
01804 #ifdef PETSC_USE_64BIT_INDICES
01805 PetscInt begin, end;
01806 #else
01807 int begin, end;
01808 #endif
01809 const int ierr = MatGetOwnershipRange (static_cast<const Mat &>(matrix),
01810 &begin, &end);
01811 AssertThrow (ierr == 0, ExcPETScError(ierr));
01812
01813 return ((index >= static_cast<unsigned int>(begin)) &&
01814 (index < static_cast<unsigned int>(end)));
01815 }
01816
01817
01818
01819 inline
01820 void
01821 MatrixBase::prepare_action(const LastAction::Values new_action)
01822 {
01823 if (last_action == new_action)
01824 ;
01825 else if (last_action == LastAction::none)
01826 last_action = new_action;
01827 else
01828 Assert (false, ExcWrongMode (last_action, new_action));
01829 }
01830
01831
01832
01833 inline
01834 void
01835 MatrixBase::prepare_add()
01836 {
01837 prepare_action(LastAction::add);
01838 }
01839
01840
01841
01842 inline
01843 void
01844 MatrixBase::prepare_set()
01845 {
01846 prepare_action(LastAction::insert);
01847 }
01848
01849 #endif // DOXYGEN
01850 }
01851
01852
01853 DEAL_II_NAMESPACE_CLOSE
01854
01855
01856 #endif // DEAL_II_USE_PETSC
01857
01858
01859
01860
01861 #endif
01862
01863