include/deal.II/lac/sparse_matrix.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: sparse_matrix.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 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__sparse_matrix_h
00013 #define __deal2__sparse_matrix_h
00014 
00015 
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/base/subscriptor.h>
00018 #include <deal.II/base/smartpointer.h>
00019 #include <deal.II/lac/sparsity_pattern.h>
00020 #include <deal.II/lac/identity_matrix.h>
00021 #include <deal.II/lac/exceptions.h>
00022 
00023 DEAL_II_NAMESPACE_OPEN
00024 
00025 template <typename number> class Vector;
00026 template <typename number> class FullMatrix;
00027 template <typename Matrix> class BlockMatrixBase;
00028 template <typename number> class SparseILU;
00029 
00035 namespace SparseMatrixIterators
00036 {
00037                                    // forward declaration
00038   template <typename number, bool Constness>
00039   class Iterator;
00040 
00058   template <typename number, bool Constness>
00059   class Accessor : public SparsityPatternIterators::Accessor
00060   {
00061     public:
00065       number value() const;
00066 
00070       number& value();
00071 
00079       const SparseMatrix<number>& get_matrix () const;
00080   };
00081 
00082 
00083 
00094   template <typename number>
00095   class Accessor<number,true> : public SparsityPatternIterators::Accessor
00096   {
00097     public:
00103       typedef const SparseMatrix<number> MatrixType;
00104 
00108       Accessor (MatrixType         *matrix,
00109                 const unsigned int  row,
00110                 const unsigned int  index);
00111 
00116       Accessor (MatrixType         *matrix);
00117 
00123       Accessor (const SparseMatrixIterators::Accessor<number,false> &a);
00124 
00128       number value() const;
00129 
00137       MatrixType & get_matrix () const;
00138 
00139     private:
00143       MatrixType *matrix;
00144 
00149       using SparsityPatternIterators::Accessor::advance;
00150 
00154       template <typename, bool>
00155       friend class Iterator;
00156   };
00157 
00158 
00169   template <typename number>
00170   class Accessor<number,false> : public SparsityPatternIterators::Accessor
00171   {
00172     private:
00216       class Reference
00217       {
00218         public:
00224           Reference (const Accessor *accessor,
00225                      const bool dummy);
00226 
00231           operator number () const;
00232 
00237           const Reference & operator = (const number n) const;
00238 
00243           const Reference & operator += (const number n) const;
00244 
00250           const Reference & operator -= (const number n) const;
00251 
00257           const Reference & operator *= (const number n) const;
00258 
00264           const Reference & operator /= (const number n) const;
00265 
00266         private:
00272           const Accessor *accessor;
00273       };
00274 
00275     public:
00281       typedef SparseMatrix<number> MatrixType;
00282 
00286       Accessor (MatrixType         *matrix,
00287                 const unsigned int  row,
00288                 const unsigned int  index);
00289 
00294       Accessor (MatrixType         *matrix);
00295 
00301       Reference value() const;
00302 
00310       MatrixType & get_matrix () const;
00311 
00312     private:
00316       MatrixType *matrix;
00317 
00322       using SparsityPatternIterators::Accessor::advance;
00323 
00327       template <typename, bool>
00328       friend class Iterator;
00329 
00335 #ifdef DEAL_II_NESTED_CLASS_FRIEND_BUG
00336       friend class Reference;
00337 #endif
00338   };
00339 
00340 
00341 
00356   template <typename number, bool Constness>
00357   class Iterator
00358   {
00359     public:
00365       typedef
00366       typename Accessor<number,Constness>::MatrixType
00367       MatrixType;
00368 
00374       Iterator (MatrixType        *matrix,
00375                 const unsigned int row,
00376                 const unsigned int index);
00377 
00382       Iterator (MatrixType *matrix);
00383 
00389       Iterator (const SparseMatrixIterators::Iterator<number,false> &i);
00390 
00394       Iterator & operator++ ();
00395 
00399       Iterator operator++ (int);
00400 
00404       const Accessor<number,Constness> & operator* () const;
00405 
00409       const Accessor<number,Constness> * operator-> () const;
00410 
00417       bool operator == (const Iterator &) const;
00418 
00422       bool operator != (const Iterator &) const;
00423 
00435       bool operator < (const Iterator &) const;
00436 
00442       bool operator > (const Iterator &) const;
00443 
00444     private:
00449       Accessor<number,Constness> accessor;
00450   };
00451 
00452 }
00453 
00454 
00455 //TODO: Add multithreading to the other vmult functions.
00456 
00473 template <typename number>
00474 class SparseMatrix : public virtual Subscriptor
00475 {
00476   public:
00481     typedef number value_type;
00482 
00502     typedef typename numbers::NumberTraits<number>::real_type real_type;
00503 
00511     typedef
00512     SparseMatrixIterators::Iterator<number,true>
00513     const_iterator;
00514 
00524     typedef
00525     SparseMatrixIterators::Iterator<number,false>
00526     iterator;
00527 
00538     struct Traits
00539     {
00545         static const bool zero_addition_can_be_elided = true;
00546     };
00547 
00567     SparseMatrix ();
00568 
00581     SparseMatrix (const SparseMatrix &);
00582 
00609     explicit SparseMatrix (const SparsityPattern &sparsity);
00610 
00623     SparseMatrix (const SparsityPattern &sparsity,
00624                   const IdentityMatrix  &id);
00625 
00631     virtual ~SparseMatrix ();
00632 
00652     SparseMatrix<number>& operator = (const SparseMatrix<number> &);
00653 
00666     SparseMatrix<number> &
00667     operator= (const IdentityMatrix  &id);
00668 
00685     SparseMatrix & operator = (const double d);
00686 
00711     virtual void reinit (const SparsityPattern &sparsity);
00712 
00721     virtual void clear ();
00723 
00734     bool empty () const;
00735 
00742     unsigned int m () const;
00743 
00750     unsigned int n () const;
00751 
00756     unsigned int get_row_length (const unsigned int row) const;
00757 
00767     unsigned int n_nonzero_elements () const;
00768 
00784     unsigned int n_actually_nonzero_elements (const double threshold = 0.) const;
00785 
00798     const SparsityPattern & get_sparsity_pattern () const;
00799 
00806     std::size_t memory_consumption () const;
00807 
00809 
00822     void set (const unsigned int i,
00823               const unsigned int j,
00824               const number value);
00825 
00853     template <typename number2>
00854     void set (const std::vector<unsigned int> &indices,
00855               const FullMatrix<number2>       &full_matrix,
00856               const bool                       elide_zero_values = false);
00857 
00865     template <typename number2>
00866     void set (const std::vector<unsigned int> &row_indices,
00867               const std::vector<unsigned int> &col_indices,
00868               const FullMatrix<number2>       &full_matrix,
00869               const bool                       elide_zero_values = false);
00870 
00889     template <typename number2>
00890     void set (const unsigned int               row,
00891               const std::vector<unsigned int> &col_indices,
00892               const std::vector<number2>      &values,
00893               const bool                       elide_zero_values = false);
00894 
00911     template <typename number2>
00912     void set (const unsigned int  row,
00913               const unsigned int  n_cols,
00914               const unsigned int *col_indices,
00915               const number2      *values,
00916               const bool          elide_zero_values = false);
00917 
00927     void add (const unsigned int i,
00928               const unsigned int j,
00929               const number value);
00930 
00958     template <typename number2>
00959     void add (const std::vector<unsigned int> &indices,
00960               const FullMatrix<number2>       &full_matrix,
00961               const bool                       elide_zero_values = true);
00962 
00970     template <typename number2>
00971     void add (const std::vector<unsigned int> &row_indices,
00972               const std::vector<unsigned int> &col_indices,
00973               const FullMatrix<number2>       &full_matrix,
00974               const bool                       elide_zero_values = true);
00975 
00993     template <typename number2>
00994     void add (const unsigned int               row,
00995               const std::vector<unsigned int> &col_indices,
00996               const std::vector<number2>      &values,
00997               const bool                       elide_zero_values = true);
00998 
01016     template <typename number2>
01017     void add (const unsigned int  row,
01018               const unsigned int  n_cols,
01019               const unsigned int *col_indices,
01020               const number2      *values,
01021               const bool          elide_zero_values = true,
01022               const bool          col_indices_are_sorted = false);
01023 
01028     SparseMatrix & operator *= (const number factor);
01029 
01034     SparseMatrix & operator /= (const number factor);
01035 
01061     void symmetrize ();
01062 
01089     template <typename somenumber>
01090     SparseMatrix<number> &
01091     copy_from (const SparseMatrix<somenumber> &source);
01092 
01127     template <typename ForwardIterator>
01128     void copy_from (const ForwardIterator begin,
01129                     const ForwardIterator end);
01130 
01141     template <typename somenumber>
01142     void copy_from (const FullMatrix<somenumber> &matrix);
01143 
01161     template <typename somenumber>
01162     void add (const number factor,
01163               const SparseMatrix<somenumber> &matrix);
01164 
01166 
01170 
01195     number operator () (const unsigned int i,
01196                         const unsigned int j) const;
01197 
01222     number el (const unsigned int i,
01223                const unsigned int j) const;
01224 
01242     number diag_element (const unsigned int i) const;
01243 
01249     number & diag_element (const unsigned int i);
01250 
01265     number raw_entry (const unsigned int row,
01266                       const unsigned int index) const;
01267 
01291     number global_entry (const unsigned int i) const;
01292 
01301     number & global_entry (const unsigned int i);
01302 
01304 
01329     template <class OutVector, class InVector>
01330     void vmult (OutVector& dst,
01331                 const InVector& src) const;
01332 
01357     template <class OutVector, class InVector>
01358     void Tvmult (OutVector& dst,
01359                  const InVector& src) const;
01360 
01384     template <class OutVector, class InVector>
01385     void vmult_add (OutVector& dst,
01386                     const InVector& src) const;
01387 
01414     template <class OutVector, class InVector>
01415     void Tvmult_add (OutVector& dst,
01416                      const InVector& src) const;
01417 
01445     template <typename somenumber>
01446     somenumber matrix_norm_square (const Vector<somenumber> &v) const;
01447 
01452     template <typename somenumber>
01453     somenumber matrix_scalar_product (const Vector<somenumber> &u,
01454                                       const Vector<somenumber> &v) const;
01455 
01471     template <typename somenumber>
01472     somenumber residual (Vector<somenumber>       &dst,
01473                          const Vector<somenumber> &x,
01474                          const Vector<somenumber> &b) const;
01475 
01516     template <typename numberB, typename numberC>
01517     void mmult (SparseMatrix<numberC>       &C,
01518                 const SparseMatrix<numberB> &B,
01519                 const Vector<number>        &V = Vector<number>(),
01520                 const bool                   rebuild_sparsity_pattern = true) const;
01521 
01564     template <typename numberB, typename numberC>
01565     void Tmmult (SparseMatrix<numberC>       &C,
01566                  const SparseMatrix<numberB> &B,
01567                  const Vector<number>       &V = Vector<number>(),
01568                  const bool                   rebuild_sparsity_pattern = true) const;
01569 
01571 
01575 
01588     real_type l1_norm () const;
01589 
01603     real_type linfty_norm () const;
01604 
01611     real_type frobenius_norm () const;
01613 
01617 
01628     template <typename somenumber>
01629     void precondition_Jacobi (Vector<somenumber>       &dst,
01630                               const Vector<somenumber> &src,
01631                               const number              omega = 1.) const;
01632 
01644     template <typename somenumber>
01645     void precondition_SSOR (Vector<somenumber>             &dst,
01646                             const Vector<somenumber>       &src,
01647                             const number                    omega = 1.,
01648                             const std::vector<unsigned int>&pos_right_of_diagonal=std::vector<unsigned int>()) const;
01649 
01654     template <typename somenumber>
01655     void precondition_SOR (Vector<somenumber>       &dst,
01656                            const Vector<somenumber> &src,
01657                            const number              om = 1.) const;
01658 
01664     template <typename somenumber>
01665     void precondition_TSOR (Vector<somenumber>       &dst,
01666                             const Vector<somenumber> &src,
01667                             const number              om = 1.) const;
01668 
01677     template <typename somenumber>
01678     void SSOR (Vector<somenumber> &v,
01679                const number        omega = 1.) const;
01680 
01686     template <typename somenumber>
01687     void SOR (Vector<somenumber> &v,
01688               const number        om = 1.) const;
01689 
01696     template <typename somenumber>
01697     void TSOR (Vector<somenumber> &v,
01698               const number        om = 1.) const;
01699 
01717     template <typename somenumber>
01718     void PSOR (Vector<somenumber> &v,
01719               const std::vector<unsigned int>& permutation,
01720               const std::vector<unsigned int>& inverse_permutation,
01721               const number        om = 1.) const;
01722 
01743     template <typename somenumber>
01744     void TPSOR (Vector<somenumber> &v,
01745               const std::vector<unsigned int>& permutation,
01746               const std::vector<unsigned int>& inverse_permutation,
01747               const number        om = 1.) const;
01748 
01758     template <typename somenumber>
01759     void Jacobi_step (Vector<somenumber> &v,
01760                       const Vector<somenumber> &b,
01761                       const number        om = 1.) const;
01762 
01769     template <typename somenumber>
01770     void SOR_step (Vector<somenumber> &v,
01771                    const Vector<somenumber> &b,
01772                    const number        om = 1.) const;
01773 
01780     template <typename somenumber>
01781     void TSOR_step (Vector<somenumber> &v,
01782                     const Vector<somenumber> &b,
01783                     const number        om = 1.) const;
01784 
01792     template <typename somenumber>
01793     void SSOR_step (Vector<somenumber> &v,
01794                     const Vector<somenumber> &b,
01795                     const number        om = 1.) const;
01797 
01801 
01807     const_iterator begin () const;
01808 
01813     const_iterator end () const;
01814 
01820     iterator begin ();
01821 
01826     iterator end ();
01827 
01841     const_iterator begin (const unsigned int r) const;
01842 
01855     const_iterator end (const unsigned int r) const;
01856 
01870     iterator begin (const unsigned int r);
01871 
01884     iterator end (const unsigned int r);
01886 
01890 
01901     template <class STREAM>
01902     void print (STREAM &out, bool across=false) const;
01903 
01945     void print_formatted (std::ostream       &out,
01946                           const unsigned int  precision   = 3,
01947                           const bool          scientific  = true,
01948                           const unsigned int  width       = 0,
01949                           const char         *zero_string = " ",
01950                           const double        denominator = 1.) const;
01951 
01961     void print_pattern(std::ostream& out,
01962                        const double threshold = 0.) const;
01963 
01984     void block_write (std::ostream &out) const;
01985 
02017     void block_read (std::istream &in);
02019 
02025     DeclException2 (ExcInvalidIndex,
02026                     int, int,
02027                     << "The entry with index <" << arg1 << ',' << arg2
02028                     << "> does not exist.");
02032     DeclException1 (ExcInvalidIndex1,
02033                     int,
02034                     << "The index " << arg1 << " is not in the allowed range.");
02038     DeclException0 (ExcDifferentSparsityPatterns);
02042     DeclException2 (ExcIteratorRange,
02043                     int, int,
02044                     << "The iterators denote a range of " << arg1
02045                     << " elements, but the given number of rows was " << arg2);
02049     DeclException0 (ExcSourceEqualsDestination);
02051 
02052   protected:
02075     void prepare_add();
02076 
02084     void prepare_set();
02085 
02086   private:
02096     SmartPointer<const SparsityPattern,SparseMatrix<number> > cols;
02097 
02110     number *val;
02111 
02123     std::size_t max_len;
02124 
02125                                      // make all other sparse matrices
02126                                      // friends
02127     template <typename somenumber> friend class SparseMatrix;
02128     template <typename somenumber> friend class SparseLUDecomposition;
02129     template <typename> friend class SparseILU;
02130 
02136     template <typename> friend class BlockMatrixBase;
02137 };
02138 
02141 #ifndef DOXYGEN
02142 /*---------------------- Inline functions -----------------------------------*/
02143 
02144 
02145 
02146 template <typename number>
02147 inline
02148 unsigned int SparseMatrix<number>::m () const
02149 {
02150   Assert (cols != 0, ExcNotInitialized());
02151   return cols->rows;
02152 }
02153 
02154 
02155 template <typename number>
02156 inline
02157 unsigned int SparseMatrix<number>::n () const
02158 {
02159   Assert (cols != 0, ExcNotInitialized());
02160   return cols->cols;
02161 }
02162 
02163 
02164                                         // Inline the set() and add()
02165                                         // functions, since they will be
02166                                         // called frequently.
02167 template <typename number>
02168 inline
02169 void
02170 SparseMatrix<number>::set (const unsigned int i,
02171                            const unsigned int j,
02172                            const number       value)
02173 {
02174   Assert (numbers::is_finite(value), ExcNumberNotFinite());
02175 
02176   const unsigned int index = cols->operator()(i, j);
02177 
02178                                    // it is allowed to set elements of
02179                                    // the matrix that are not part of
02180                                    // the sparsity pattern, if the
02181                                    // value to which we set it is zero
02182   if (index == SparsityPattern::invalid_entry)
02183     {
02184       Assert ((index != SparsityPattern::invalid_entry) ||
02185               (value == 0.),
02186               ExcInvalidIndex(i, j));
02187       return;
02188     }
02189 
02190   val[index] = value;
02191 }
02192 
02193 
02194 
02195 template <typename number>
02196 template <typename number2>
02197 inline
02198 void
02199 SparseMatrix<number>::set (const std::vector<unsigned int> &indices,
02200                            const FullMatrix<number2>       &values,
02201                            const bool                       elide_zero_values)
02202 {
02203   Assert (indices.size() == values.m(),
02204           ExcDimensionMismatch(indices.size(), values.m()));
02205   Assert (values.m() == values.n(), ExcNotQuadratic());
02206 
02207   for (unsigned int i=0; i<indices.size(); ++i)
02208     set (indices[i], indices.size(), &indices[0], &values(i,0),
02209          elide_zero_values);
02210 }
02211 
02212 
02213 
02214 template <typename number>
02215 template <typename number2>
02216 inline
02217 void
02218 SparseMatrix<number>::set (const std::vector<unsigned int> &row_indices,
02219                            const std::vector<unsigned int> &col_indices,
02220                            const FullMatrix<number2>       &values,
02221                            const bool                       elide_zero_values)
02222 {
02223   Assert (row_indices.size() == values.m(),
02224           ExcDimensionMismatch(row_indices.size(), values.m()));
02225   Assert (col_indices.size() == values.n(),
02226           ExcDimensionMismatch(col_indices.size(), values.n()));
02227 
02228   for (unsigned int i=0; i<row_indices.size(); ++i)
02229     set (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
02230          elide_zero_values);
02231 }
02232 
02233 
02234 
02235 template <typename number>
02236 template <typename number2>
02237 inline
02238 void
02239 SparseMatrix<number>::set (const unsigned int               row,
02240                            const std::vector<unsigned int> &col_indices,
02241                            const std::vector<number2>      &values,
02242                            const bool                       elide_zero_values)
02243 {
02244   Assert (col_indices.size() == values.size(),
02245           ExcDimensionMismatch(col_indices.size(), values.size()));
02246 
02247   set (row, col_indices.size(), &col_indices[0], &values[0],
02248        elide_zero_values);
02249 }
02250 
02251 
02252 
02253 template <typename number>
02254 inline
02255 void
02256 SparseMatrix<number>::add (const unsigned int i,
02257                            const unsigned int j,
02258                            const number       value)
02259 {
02260   Assert (numbers::is_finite(value), ExcNumberNotFinite());
02261 
02262   if (value == 0)
02263     return;
02264 
02265   const unsigned int index = cols->operator()(i, j);
02266 
02267                                    // it is allowed to add elements to
02268                                    // the matrix that are not part of
02269                                    // the sparsity pattern, if the
02270                                    // value to which we set it is zero
02271   if (index == SparsityPattern::invalid_entry)
02272     {
02273       Assert ((index != SparsityPattern::invalid_entry) ||
02274               (value == 0.),
02275               ExcInvalidIndex(i, j));
02276       return;
02277     }
02278 
02279   val[index] += value;
02280 }
02281 
02282 
02283 
02284 template <typename number>
02285 template <typename number2>
02286 inline
02287 void
02288 SparseMatrix<number>::add (const std::vector<unsigned int> &indices,
02289                            const FullMatrix<number2>       &values,
02290                            const bool                       elide_zero_values)
02291 {
02292   Assert (indices.size() == values.m(),
02293           ExcDimensionMismatch(indices.size(), values.m()));
02294   Assert (values.m() == values.n(), ExcNotQuadratic());
02295 
02296   for (unsigned int i=0; i<indices.size(); ++i)
02297     add (indices[i], indices.size(), &indices[0], &values(i,0),
02298          elide_zero_values);
02299 }
02300 
02301 
02302 
02303 template <typename number>
02304 template <typename number2>
02305 inline
02306 void
02307 SparseMatrix<number>::add (const std::vector<unsigned int> &row_indices,
02308                            const std::vector<unsigned int> &col_indices,
02309                            const FullMatrix<number2>       &values,
02310                            const bool                       elide_zero_values)
02311 {
02312   Assert (row_indices.size() == values.m(),
02313           ExcDimensionMismatch(row_indices.size(), values.m()));
02314   Assert (col_indices.size() == values.n(),
02315           ExcDimensionMismatch(col_indices.size(), values.n()));
02316 
02317   for (unsigned int i=0; i<row_indices.size(); ++i)
02318     add (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
02319          elide_zero_values);
02320 }
02321 
02322 
02323 
02324 template <typename number>
02325 template <typename number2>
02326 inline
02327 void
02328 SparseMatrix<number>::add (const unsigned int               row,
02329                            const std::vector<unsigned int> &col_indices,
02330                            const std::vector<number2>      &values,
02331                            const bool                       elide_zero_values)
02332 {
02333   Assert (col_indices.size() == values.size(),
02334           ExcDimensionMismatch(col_indices.size(), values.size()));
02335 
02336   add (row, col_indices.size(), &col_indices[0], &values[0],
02337        elide_zero_values);
02338 }
02339 
02340 
02341 
02342 template <typename number>
02343 inline
02344 SparseMatrix<number> &
02345 SparseMatrix<number>::operator *= (const number factor)
02346 {
02347   Assert (cols != 0, ExcNotInitialized());
02348   Assert (val != 0, ExcNotInitialized());
02349 
02350   number             *val_ptr    = &val[0];
02351   const number *const end_ptr    = &val[cols->n_nonzero_elements()];
02352 
02353   while (val_ptr != end_ptr)
02354     *val_ptr++ *= factor;
02355 
02356   return *this;
02357 }
02358 
02359 
02360 
02361 template <typename number>
02362 inline
02363 SparseMatrix<number> &
02364 SparseMatrix<number>::operator /= (const number factor)
02365 {
02366   Assert (cols != 0, ExcNotInitialized());
02367   Assert (val != 0, ExcNotInitialized());
02368   Assert (factor !=0, ExcDivideByZero());
02369 
02370   const number factor_inv = 1. / factor;
02371 
02372   number             *val_ptr    = &val[0];
02373   const number *const end_ptr    = &val[cols->n_nonzero_elements()];
02374 
02375   while (val_ptr != end_ptr)
02376     *val_ptr++ *= factor_inv;
02377 
02378   return *this;
02379 }
02380 
02381 
02382 
02383 template <typename number>
02384 inline
02385 number SparseMatrix<number>::operator () (const unsigned int i,
02386                                           const unsigned int j) const
02387 {
02388   Assert (cols != 0, ExcNotInitialized());
02389   Assert (cols->operator()(i,j) != SparsityPattern::invalid_entry,
02390                ExcInvalidIndex(i,j));
02391   return val[cols->operator()(i,j)];
02392 }
02393 
02394 
02395 
02396 template <typename number>
02397 inline
02398 number SparseMatrix<number>::el (const unsigned int i,
02399                                  const unsigned int j) const
02400 {
02401   Assert (cols != 0, ExcNotInitialized());
02402   const unsigned int index = cols->operator()(i,j);
02403 
02404   if (index != SparsityPattern::invalid_entry)
02405     return val[index];
02406   else
02407     return 0;
02408 }
02409 
02410 
02411 
02412 template <typename number>
02413 inline
02414 number SparseMatrix<number>::diag_element (const unsigned int i) const
02415 {
02416   Assert (cols != 0, ExcNotInitialized());
02417   Assert (cols->optimize_diagonal(),  ExcNotQuadratic());
02418   Assert (i<m(), ExcInvalidIndex1(i));
02419 
02420                                    // Use that the first element in each
02421                                    // row of a quadratic matrix is the main
02422                                    // diagonal
02423   return val[cols->rowstart[i]];
02424 }
02425 
02426 
02427 
02428 template <typename number>
02429 inline
02430 number & SparseMatrix<number>::diag_element (const unsigned int i)
02431 {
02432   Assert (cols != 0, ExcNotInitialized());
02433   Assert (cols->optimize_diagonal(),  ExcNotQuadratic());
02434   Assert (i<m(), ExcInvalidIndex1(i));
02435 
02436                                    // Use that the first element in each
02437                                    // row of a quadratic matrix is the main
02438                                    // diagonal
02439   return val[cols->rowstart[i]];
02440 }
02441 
02442 
02443 
02444 template <typename number>
02445 inline
02446 number
02447 SparseMatrix<number>::raw_entry (const unsigned int row,
02448                                  const unsigned int index) const
02449 {
02450   Assert(row<cols->rows, ExcIndexRange(row,0,cols->rows));
02451   Assert(index<cols->row_length(row),
02452          ExcIndexRange(index,0,cols->row_length(row)));
02453 
02454   return val[cols->rowstart[row]+index];
02455 }
02456 
02457 
02458 
02459 template <typename number>
02460 inline
02461 number SparseMatrix<number>::global_entry (const unsigned int j) const
02462 {
02463   Assert (cols != 0, ExcNotInitialized());
02464   Assert (j < cols->n_nonzero_elements(),
02465           ExcIndexRange (j, 0, cols->n_nonzero_elements()));
02466 
02467   return val[j];
02468 }
02469 
02470 
02471 
02472 template <typename number>
02473 inline
02474 number & SparseMatrix<number>::global_entry (const unsigned int j)
02475 {
02476   Assert (cols != 0, ExcNotInitialized());
02477   Assert (j < cols->n_nonzero_elements(),
02478           ExcIndexRange (j, 0, cols->n_nonzero_elements()));
02479 
02480   return val[j];
02481 }
02482 
02483 
02484 
02485 template <typename number>
02486 template <typename ForwardIterator>
02487 void
02488 SparseMatrix<number>::copy_from (const ForwardIterator begin,
02489                                  const ForwardIterator end)
02490 {
02491   Assert (static_cast<unsigned int>(std::distance (begin, end)) == m(),
02492           ExcIteratorRange (std::distance (begin, end), m()));
02493 
02494                                    // for use in the inner loop, we
02495                                    // define a typedef to the type of
02496                                    // the inner iterators
02497   typedef typename std::iterator_traits<ForwardIterator>::value_type::const_iterator inner_iterator;
02498   unsigned int row=0;
02499   for (ForwardIterator i=begin; i!=end; ++i, ++row)
02500     {
02501       const inner_iterator end_of_row = i->end();
02502       for (inner_iterator j=i->begin(); j!=end_of_row; ++j)
02503                                          // write entries
02504         set (row, j->first, j->second);
02505     };
02506 }
02507 
02508 
02509 //---------------------------------------------------------------------------
02510 
02511 
02512 namespace SparseMatrixIterators
02513 {
02514   template <typename number>
02515   inline
02516   Accessor<number,true>::
02517   Accessor (const MatrixType   *matrix,
02518             const unsigned int  row,
02519             const unsigned int  index)
02520                   :
02521                   SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
02522                                                       row, index),
02523                   matrix (matrix)
02524   {}
02525 
02526 
02527 
02528   template <typename number>
02529   inline
02530   Accessor<number,true>::
02531   Accessor (const MatrixType *matrix)
02532                   :
02533                   SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern()),
02534                   matrix (matrix)
02535   {}
02536 
02537 
02538 
02539   template <typename number>
02540   inline
02541   Accessor<number,true>::
02542   Accessor (const SparseMatrixIterators::Accessor<number,false> &a)
02543                   :
02544                   SparsityPatternIterators::Accessor (a),
02545                   matrix (&a.get_matrix())
02546   {}
02547 
02548 
02549 
02550   template <typename number>
02551   inline
02552   number
02553   Accessor<number, true>::value () const
02554   {
02555     return matrix->raw_entry(a_row, a_index);
02556   }
02557 
02558 
02559 
02560   template <typename number>
02561   inline
02562   typename Accessor<number, true>::MatrixType &
02563   Accessor<number, true>::get_matrix () const
02564   {
02565     return *matrix;
02566   }
02567 
02568 
02569 
02570   template <typename number>
02571   inline
02572   Accessor<number, false>::Reference::
02573   Reference (const Accessor *accessor,
02574              const bool)
02575                   :
02576                   accessor (accessor)
02577   {}
02578 
02579 
02580 
02581   template <typename number>
02582   inline
02583   Accessor<number, false>::Reference::operator number() const
02584   {
02585     return accessor->matrix->raw_entry(accessor->a_row,
02586                                        accessor->a_index);
02587   }
02588 
02589 
02590 
02591   template <typename number>
02592   inline
02593   const typename Accessor<number, false>::Reference &
02594   Accessor<number, false>::Reference::operator = (const number n) const
02595   {
02596 //TODO: one could optimize this by not going again through the mapping from row/col index to global index
02597     accessor->matrix->set (accessor->row(), accessor->column(), n);
02598     return *this;
02599   }
02600 
02601 
02602 
02603   template <typename number>
02604   inline
02605   const typename Accessor<number, false>::Reference &
02606   Accessor<number, false>::Reference::operator += (const number n) const
02607   {
02608 //TODO: one could optimize this by not going again through the mapping from row/col index to global index
02609     accessor->matrix->set (accessor->row(), accessor->column(),
02610                            static_cast<number>(*this) + n);
02611     return *this;
02612   }
02613 
02614 
02615 
02616   template <typename number>
02617   inline
02618   const typename Accessor<number, false>::Reference &
02619   Accessor<number, false>::Reference::operator -= (const number n) const
02620   {
02621 //TODO: one could optimize this by not going again through the mapping from row/col index to global index
02622     accessor->matrix->set (accessor->row(), accessor->column(),
02623                            static_cast<number>(*this) - n);
02624     return *this;
02625   }
02626 
02627 
02628 
02629   template <typename number>
02630   inline
02631   const typename Accessor<number, false>::Reference &
02632   Accessor<number, false>::Reference::operator *= (const number n) const
02633   {
02634 //TODO: one could optimize this by not going again through the mapping from row/col index to global index
02635     accessor->matrix->set (accessor->row(), accessor->column(),
02636                            static_cast<number>(*this)*n);
02637     return *this;
02638   }
02639 
02640 
02641 
02642   template <typename number>
02643   inline
02644   const typename Accessor<number, false>::Reference &
02645   Accessor<number, false>::Reference::operator /= (const number n) const
02646   {
02647 //TODO: one could optimize this by not going again through the mapping from row/col index to global index
02648     accessor->matrix->set (accessor->row(), accessor->column(),
02649                            static_cast<number>(*this)/n);
02650     return *this;
02651   }
02652 
02653 
02654 
02655   template <typename number>
02656   inline
02657   Accessor<number,false>::
02658   Accessor (MatrixType         *matrix,
02659             const unsigned int  row,
02660             const unsigned int  index)
02661                   :
02662                   SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
02663                                                       row, index),
02664                   matrix (matrix)
02665   {}
02666 
02667 
02668 
02669   template <typename number>
02670   inline
02671   Accessor<number,false>::
02672   Accessor (MatrixType         *matrix)
02673                   :
02674                   SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern()),
02675                   matrix (matrix)
02676   {}
02677 
02678 
02679 
02680   template <typename number>
02681   inline
02682   typename Accessor<number, false>::Reference
02683   Accessor<number, false>::value() const
02684   {
02685     return Reference(this,true);
02686   }
02687 
02688 
02689 
02690 
02691   template <typename number>
02692   inline
02693   typename Accessor<number, false>::MatrixType &
02694   Accessor<number, false>::get_matrix () const
02695   {
02696     return *matrix;
02697   }
02698 
02699 
02700 
02701   template <typename number, bool Constness>
02702   inline
02703   Iterator<number, Constness>::
02704   Iterator (MatrixType        *matrix,
02705             const unsigned int r,
02706             const unsigned int i)
02707                   :
02708                   accessor(matrix, r, i)
02709   {}
02710 
02711 
02712 
02713   template <typename number, bool Constness>
02714   inline
02715   Iterator<number, Constness>::
02716   Iterator (MatrixType *matrix)
02717                   :
02718                   accessor(matrix)
02719   {}
02720 
02721 
02722 
02723   template <typename number, bool Constness>
02724   inline
02725   Iterator<number, Constness>::
02726   Iterator (const SparseMatrixIterators::Iterator<number,false> &i)
02727                   :
02728                   accessor(*i)
02729   {}
02730 
02731 
02732 
02733   template <typename number, bool Constness>
02734   inline
02735   Iterator<number, Constness> &
02736   Iterator<number,Constness>::operator++ ()
02737   {
02738     accessor.advance ();
02739     return *this;
02740   }
02741 
02742 
02743   template <typename number, bool Constness>
02744   inline
02745   Iterator<number,Constness>
02746   Iterator<number,Constness>::operator++ (int)
02747   {
02748     const Iterator iter = *this;
02749     accessor.advance ();
02750     return iter;
02751   }
02752 
02753 
02754   template <typename number, bool Constness>
02755   inline
02756   const Accessor<number,Constness> &
02757   Iterator<number,Constness>::operator* () const
02758   {
02759     return accessor;
02760   }
02761 
02762 
02763   template <typename number, bool Constness>
02764   inline
02765   const Accessor<number,Constness> *
02766   Iterator<number,Constness>::operator-> () const
02767   {
02768     return &accessor;
02769   }
02770 
02771 
02772   template <typename number, bool Constness>
02773   inline
02774   bool
02775   Iterator<number,Constness>::
02776   operator == (const Iterator& other) const
02777   {
02778     return (accessor == other.accessor);
02779   }
02780 
02781 
02782   template <typename number, bool Constness>
02783   inline
02784   bool
02785   Iterator<number,Constness>::
02786   operator != (const Iterator& other) const
02787   {
02788     return ! (*this == other);
02789   }
02790 
02791 
02792   template <typename number, bool Constness>
02793   inline
02794   bool
02795   Iterator<number,Constness>::
02796   operator < (const Iterator& other) const
02797   {
02798     Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
02799             ExcInternalError());
02800 
02801     return (accessor < other.accessor);
02802   }
02803 
02804 
02805   template <typename number, bool Constness>
02806   inline
02807   bool
02808   Iterator<number,Constness>::
02809   operator > (const Iterator& other) const
02810   {
02811     return (other < *this);
02812   }
02813 
02814 }
02815 
02816 
02817 
02818 template <typename number>
02819 inline
02820 typename SparseMatrix<number>::const_iterator
02821 SparseMatrix<number>::begin () const
02822 {
02823                                    // search for the first line with a nonzero
02824                                    // number of entries
02825   for (unsigned int r=0; r<m(); ++r)
02826     if (cols->row_length(r) > 0)
02827       return const_iterator(this, r, 0);
02828 
02829                                    // alright, this matrix is completely
02830                                    // empty. that's strange but ok. simply
02831                                    // return the end() iterator
02832   return end();
02833 }
02834 
02835 
02836 template <typename number>
02837 inline
02838 typename SparseMatrix<number>::const_iterator
02839 SparseMatrix<number>::end () const
02840 {
02841   return const_iterator(this);
02842 }
02843 
02844 
02845 template <typename number>
02846 inline
02847 typename SparseMatrix<number>::iterator
02848 SparseMatrix<number>::begin ()
02849 {
02850                                    // search for the first line with a nonzero
02851                                    // number of entries
02852   for (unsigned int r=0; r<m(); ++r)
02853     if (cols->row_length(r) > 0)
02854       return iterator(this, r, 0);
02855 
02856                                    // alright, this matrix is completely
02857                                    // empty. that's strange but ok. simply
02858                                    // return the end() iterator
02859   return end();
02860 }
02861 
02862 
02863 template <typename number>
02864 inline
02865 typename SparseMatrix<number>::iterator
02866 SparseMatrix<number>::end ()
02867 {
02868   return iterator(this, m(), 0);
02869 }
02870 
02871 
02872 template <typename number>
02873 inline
02874 typename SparseMatrix<number>::const_iterator
02875 SparseMatrix<number>::begin (const unsigned int r) const
02876 {
02877   Assert (r<m(), ExcIndexRange(r,0,m()));
02878 
02879   if (cols->row_length(r) > 0)
02880     return const_iterator(this, r, 0);
02881   else
02882     return end (r);
02883 }
02884 
02885 
02886 
02887 template <typename number>
02888 inline
02889 typename SparseMatrix<number>::const_iterator
02890 SparseMatrix<number>::end (const unsigned int r) const
02891 {
02892   Assert (r<m(), ExcIndexRange(r,0,m()));
02893 
02894                                    // place the iterator on the first entry
02895                                    // past this line, or at the end of the
02896                                    // matrix
02897   for (unsigned int i=r+1; i<m(); ++i)
02898     if (cols->row_length(i) > 0)
02899       return const_iterator(this, i, 0);
02900 
02901                                    // if there is no such line, then take the
02902                                    // end iterator of the matrix
02903   return end();
02904 }
02905 
02906 
02907 
02908 template <typename number>
02909 inline
02910 typename SparseMatrix<number>::iterator
02911 SparseMatrix<number>::begin (const unsigned int r)
02912 {
02913   Assert (r<m(), ExcIndexRange(r,0,m()));
02914 
02915   if (cols->row_length(r) > 0)
02916     return iterator(this, r, 0);
02917   else
02918     return end (r);
02919 }
02920 
02921 
02922 
02923 template <typename number>
02924 inline
02925 typename SparseMatrix<number>::iterator
02926 SparseMatrix<number>::end (const unsigned int r)
02927 {
02928   Assert (r<m(), ExcIndexRange(r,0,m()));
02929 
02930                                    // place the iterator on the first entry
02931                                    // past this line, or at the end of the
02932                                    // matrix
02933   for (unsigned int i=r+1; i<m(); ++i)
02934     if (cols->row_length(i) > 0)
02935       return iterator(this, i, 0);
02936 
02937                                    // if there is no such line, then take the
02938                                    // end iterator of the matrix
02939   return end();
02940 }
02941 
02942 
02943 
02944 template <typename number>
02945 template <class STREAM>
02946 inline
02947 void SparseMatrix<number>::print (STREAM &out, bool across) const
02948 {
02949   Assert (cols != 0, ExcNotInitialized());
02950   Assert (val != 0, ExcNotInitialized());
02951 
02952   if (across)
02953     {
02954       for (unsigned int i=0; i<cols->rows; ++i)
02955         for (unsigned int j=cols->rowstart[i]; j<cols->rowstart[i+1]; ++j)
02956           out << ' ' << i << ',' << cols->colnums[j] << ':' << val[j];
02957       out << std::endl;
02958     }
02959   else
02960     for (unsigned int i=0; i<cols->rows; ++i)
02961       for (unsigned int j=cols->rowstart[i]; j<cols->rowstart[i+1]; ++j)
02962         out << "(" << i << "," << cols->colnums[j] << ") " << val[j] << std::endl;
02963 }
02964 
02965 
02966 template <typename number>
02967 inline
02968 void
02969 SparseMatrix<number>::
02970 prepare_add()
02971 {
02972                                    //nothing to do here
02973 }
02974 
02975 
02976 
02977 template <typename number>
02978 inline
02979 void
02980 SparseMatrix<number>::
02981 prepare_set()
02982 {
02983                                    //nothing to do here
02984 }
02985 
02986 #endif // DOXYGEN
02987 
02988 
02989 /*----------------------------   sparse_matrix.h     ---------------------------*/
02990 
02991 DEAL_II_NAMESPACE_CLOSE
02992 
02993 #endif
02994 /*----------------------------   sparse_matrix.h     ---------------------------*/
02995 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

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