include/deal.II/lac/petsc_matrix_base.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: petsc_matrix_base.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 2004, 2005, 2006, 2007, 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__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                                    // forward declarations
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                                          * Abstract PETSc object that helps view
01138                                          * in ASCII other PETSc objects. Currently
01139                                          * this function simply writes non-zero
01140                                          * elements of a matrix to the terminal.
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 // -------------------------- inline and template functions ----------------------
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                                        // if at end of line: do one step, then
01374                                        // cycle until we find a row with a
01375                                        // nonzero number of entries
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                                         // Inline the set() and add()
01452                                         // functions, since they will be
01453                                         // called frequently, and the
01454                                         // compiler can optimize away
01455                                         // some unnecessary loops when
01456                                         // the sizes are given at
01457                                         // compile time.
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                                    // If we don't elide zeros, the pointers
01544                                    // are already available...
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                                    // Otherwise, extract nonzero values in
01556                                    // each row and get the respective index.
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                                   // we have to do checkings on Insert/Add
01601                                   // in any case
01602                                   // to be consistent with the MPI
01603                                   // communication model (see the comments
01604                                   // in the documentation of
01605                                   // TrilinosWrappers::Vector), but we can
01606                                   // save some work if the addend is
01607                                   // zero. However, these actions are done
01608                                   // in case we pass on to the other
01609                                   // function.
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          /*col_indices_are_sorted*/)
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                                    // If we don't elide zeros, the pointers
01694                                    // are already available...
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                                    // Otherwise, extract nonzero values in
01706                                    // each row and get the respective index.
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                                      // place the iterator on the first entry
01787                                      // past this line, or at the end of the
01788                                      // matrix
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                                      // if there is no such line, then take the
01794                                      // end iterator of the matrix
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 /*----------------------------   petsc_matrix_base.h     ---------------------------*/
01860 
01861 #endif
01862 /*----------------------------   petsc_matrix_base.h     ---------------------------*/
01863 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Wed May 23 2012 06:07:38 by doxygen 1.7.3