include/deal.II/lac/chunk_sparse_matrix.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: chunk_sparse_matrix.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 2008, 2009, 2010, 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__chunk_sparse_matrix_h
00013 #define __deal2__chunk_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/chunk_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 
00048 template <typename number>
00049 class ChunkSparseMatrix : public virtual Subscriptor
00050 {
00051   public:
00056     typedef number value_type;
00057 
00077     typedef typename numbers::NumberTraits<number>::real_type real_type;
00078 
00089     struct Traits
00090     {
00096         static const bool zero_addition_can_be_elided = true;
00097     };
00098 
00118     ChunkSparseMatrix ();
00119 
00132     ChunkSparseMatrix (const ChunkSparseMatrix &);
00133 
00160     explicit ChunkSparseMatrix (const ChunkSparsityPattern &sparsity);
00161 
00174     ChunkSparseMatrix (const ChunkSparsityPattern &sparsity,
00175                        const IdentityMatrix  &id);
00176 
00182     virtual ~ChunkSparseMatrix ();
00183 
00203     ChunkSparseMatrix<number>& operator = (const ChunkSparseMatrix<number> &);
00204 
00217     ChunkSparseMatrix<number> &
00218     operator= (const IdentityMatrix  &id);
00219 
00236     ChunkSparseMatrix & operator = (const double d);
00237 
00262     virtual void reinit (const ChunkSparsityPattern &sparsity);
00263 
00272     virtual void clear ();
00274 
00285     bool empty () const;
00286 
00293     unsigned int m () const;
00294 
00301     unsigned int n () const;
00302 
00312     unsigned int n_nonzero_elements () const;
00313 
00326     unsigned int n_actually_nonzero_elements () const;
00327 
00340     const ChunkSparsityPattern & get_sparsity_pattern () const;
00341 
00348     std::size_t memory_consumption () const;
00349 
00351 
00364     void set (const unsigned int i,
00365               const unsigned int j,
00366               const number value);
00367 
00377     void add (const unsigned int i,
00378               const unsigned int j,
00379               const number value);
00380 
00385     ChunkSparseMatrix & operator *= (const number factor);
00386 
00391     ChunkSparseMatrix & operator /= (const number factor);
00392 
00418     void symmetrize ();
00419 
00446     template <typename somenumber>
00447     ChunkSparseMatrix<number> &
00448     copy_from (const ChunkSparseMatrix<somenumber> &source);
00449 
00484     template <typename ForwardIterator>
00485     void copy_from (const ForwardIterator begin,
00486                     const ForwardIterator end);
00487 
00498     template <typename somenumber>
00499     void copy_from (const FullMatrix<somenumber> &matrix);
00500 
00518     template <typename somenumber>
00519     void add (const number factor,
00520               const ChunkSparseMatrix<somenumber> &matrix);
00521 
00523 
00527 
00552     number operator () (const unsigned int i,
00553                         const unsigned int j) const;
00554 
00579     number el (const unsigned int i,
00580                const unsigned int j) const;
00581 
00599     number diag_element (const unsigned int i) const;
00600 
00606     number & diag_element (const unsigned int i);
00607 
00609 
00634     template <class OutVector, class InVector>
00635     void vmult (OutVector& dst,
00636                 const InVector& src) const;
00637 
00662     template <class OutVector, class InVector>
00663     void Tvmult (OutVector& dst,
00664                  const InVector& src) const;
00665 
00689     template <class OutVector, class InVector>
00690     void vmult_add (OutVector& dst,
00691                     const InVector& src) const;
00692 
00719     template <class OutVector, class InVector>
00720     void Tvmult_add (OutVector& dst,
00721                      const InVector& src) const;
00722 
00750     template <typename somenumber>
00751     somenumber matrix_norm_square (const Vector<somenumber> &v) const;
00752 
00757     template <typename somenumber>
00758     somenumber matrix_scalar_product (const Vector<somenumber> &u,
00759                                       const Vector<somenumber> &v) const;
00775     template <typename somenumber>
00776     somenumber residual (Vector<somenumber>       &dst,
00777                          const Vector<somenumber> &x,
00778                          const Vector<somenumber> &b) const;
00779 
00781 
00785 
00797     real_type l1_norm () const;
00798 
00811     real_type linfty_norm () const;
00812 
00819     real_type frobenius_norm () const;
00821 
00825 
00836     template <typename somenumber>
00837     void precondition_Jacobi (Vector<somenumber>       &dst,
00838                               const Vector<somenumber> &src,
00839                               const number              omega = 1.) const;
00840 
00845     template <typename somenumber>
00846     void precondition_SSOR (Vector<somenumber>       &dst,
00847                             const Vector<somenumber> &src,
00848                             const number              om = 1.) const;
00849 
00854     template <typename somenumber>
00855     void precondition_SOR (Vector<somenumber>       &dst,
00856                            const Vector<somenumber> &src,
00857                            const number              om = 1.) const;
00858 
00864     template <typename somenumber>
00865     void precondition_TSOR (Vector<somenumber>       &dst,
00866                             const Vector<somenumber> &src,
00867                             const number              om = 1.) const;
00868 
00877     template <typename somenumber>
00878     void SSOR (Vector<somenumber> &v,
00879                const number        omega = 1.) const;
00880 
00886     template <typename somenumber>
00887     void SOR (Vector<somenumber> &v,
00888               const number        om = 1.) const;
00889 
00896     template <typename somenumber>
00897     void TSOR (Vector<somenumber> &v,
00898                const number        om = 1.) const;
00899 
00917     template <typename somenumber>
00918     void PSOR (Vector<somenumber> &v,
00919                const std::vector<unsigned int>& permutation,
00920                const std::vector<unsigned int>& inverse_permutation,
00921                const number        om = 1.) const;
00922 
00943     template <typename somenumber>
00944     void TPSOR (Vector<somenumber> &v,
00945                 const std::vector<unsigned int>& permutation,
00946                 const std::vector<unsigned int>& inverse_permutation,
00947                 const number        om = 1.) const;
00948 
00955     template <typename somenumber>
00956     void SOR_step (Vector<somenumber> &v,
00957                    const Vector<somenumber> &b,
00958                    const number        om = 1.) const;
00959 
00966     template <typename somenumber>
00967     void TSOR_step (Vector<somenumber> &v,
00968                     const Vector<somenumber> &b,
00969                     const number        om = 1.) const;
00970 
00978     template <typename somenumber>
00979     void SSOR_step (Vector<somenumber> &v,
00980                     const Vector<somenumber> &b,
00981                     const number        om = 1.) const;
00983 
00987 
00995     void print (std::ostream &out) const;
00996 
01038     void print_formatted (std::ostream       &out,
01039                           const unsigned int  precision   = 3,
01040                           const bool          scientific  = true,
01041                           const unsigned int  width       = 0,
01042                           const char         *zero_string = " ",
01043                           const double        denominator = 1.) const;
01044 
01054     void print_pattern(std::ostream& out,
01055                        const double threshold = 0.) const;
01056 
01077     void block_write (std::ostream &out) const;
01078 
01110     void block_read (std::istream &in);
01112 
01118     DeclException2 (ExcInvalidIndex,
01119                     int, int,
01120                     << "The entry with index <" << arg1 << ',' << arg2
01121                     << "> does not exist.");
01125     DeclException1 (ExcInvalidIndex1,
01126                     int,
01127                     << "The index " << arg1 << " is not in the allowed range.");
01131     DeclException0 (ExcDifferentChunkSparsityPatterns);
01135     DeclException2 (ExcIteratorRange,
01136                     int, int,
01137                     << "The iterators denote a range of " << arg1
01138                     << " elements, but the given number of rows was " << arg2);
01142     DeclException0 (ExcSourceEqualsDestination);
01144   private:
01154     SmartPointer<const ChunkSparsityPattern,ChunkSparseMatrix<number> > cols;
01155 
01168     number *val;
01169 
01181     unsigned int max_len;
01182 
01187     unsigned int compute_location (const unsigned int i,
01188                                    const unsigned int j) const;
01189 
01190                                      // make all other sparse matrices
01191                                      // friends
01192     template <typename somenumber> friend class ChunkSparseMatrix;
01193 };
01194 
01197 #ifndef DOXYGEN
01198 /*---------------------- Inline functions -----------------------------------*/
01199 
01200 
01201 
01202 template <typename number>
01203 inline
01204 unsigned int ChunkSparseMatrix<number>::m () const
01205 {
01206   Assert (cols != 0, ExcNotInitialized());
01207   return cols->rows;
01208 }
01209 
01210 
01211 template <typename number>
01212 inline
01213 unsigned int ChunkSparseMatrix<number>::n () const
01214 {
01215   Assert (cols != 0, ExcNotInitialized());
01216   return cols->cols;
01217 }
01218 
01219 
01220 
01221 template <typename number>
01222 inline
01223 unsigned int
01224 ChunkSparseMatrix<number>::compute_location (const unsigned int i,
01225                                              const unsigned int j) const
01226 {
01227   const unsigned int chunk_size = cols->get_chunk_size();
01228   const unsigned int chunk_index
01229     = cols->sparsity_pattern(i/chunk_size, j/chunk_size);
01230 
01231   if (chunk_index == ChunkSparsityPattern::invalid_entry)
01232     return ChunkSparsityPattern::invalid_entry;
01233   else
01234     {
01235       return (chunk_index * chunk_size * chunk_size
01236               +
01237               (i % chunk_size) * chunk_size
01238               +
01239               (j % chunk_size));
01240     }
01241 }
01242 
01243 
01244 template <typename number>
01245 inline
01246 void ChunkSparseMatrix<number>::set (const unsigned int i,
01247                                      const unsigned int j,
01248                                      const number value)
01249 {
01250 
01251   Assert (numbers::is_finite(value), ExcNumberNotFinite());
01252 
01253   Assert (cols != 0, ExcNotInitialized());
01254                                    // it is allowed to set elements of
01255                                    // the matrix that are not part of
01256                                    // the sparsity pattern, if the
01257                                    // value to which we set it is zero
01258   const unsigned int index = compute_location(i,j);
01259   Assert ((index != SparsityPattern::invalid_entry) ||
01260           (value == 0.),
01261           ExcInvalidIndex(i,j));
01262 
01263   if (index != SparsityPattern::invalid_entry)
01264     val[index] = value;
01265 }
01266 
01267 
01268 
01269 template <typename number>
01270 inline
01271 void ChunkSparseMatrix<number>::add (const unsigned int i,
01272                                      const unsigned int j,
01273                                      const number value)
01274 {
01275 
01276   Assert (numbers::is_finite(value), ExcNumberNotFinite());
01277 
01278   Assert (cols != 0, ExcNotInitialized());
01279 
01280   const unsigned int index = compute_location(i,j);
01281   Assert ((index != ChunkSparsityPattern::invalid_entry) ||
01282           (value == 0.),
01283           ExcInvalidIndex(i,j));
01284 
01285   if (value != 0.)
01286     val[index] += value;
01287 }
01288 
01289 
01290 
01291 template <typename number>
01292 inline
01293 ChunkSparseMatrix<number> &
01294 ChunkSparseMatrix<number>::operator *= (const number factor)
01295 {
01296   Assert (cols != 0, ExcNotInitialized());
01297   Assert (val != 0, ExcNotInitialized());
01298 
01299   const unsigned int chunk_size = cols->get_chunk_size();
01300 
01301                                    // multiply all elements of the matrix with
01302                                    // the given factor. this includes the
01303                                    // padding elements in chunks that overlap
01304                                    // the boundaries of the actual matrix --
01305                                    // but since multiplication with a number
01306                                    // does not violate the invariant of
01307                                    // keeping these elements at zero nothing
01308                                    // can happen
01309   number             *val_ptr    = val;
01310   const number *const end_ptr    = val +
01311                                    cols->sparsity_pattern.n_nonzero_elements()
01312                                    *
01313                                    chunk_size * chunk_size;
01314   while (val_ptr != end_ptr)
01315     *val_ptr++ *= factor;
01316 
01317   return *this;
01318 }
01319 
01320 
01321 
01322 template <typename number>
01323 inline
01324 ChunkSparseMatrix<number> &
01325 ChunkSparseMatrix<number>::operator /= (const number factor)
01326 {
01327   Assert (cols != 0, ExcNotInitialized());
01328   Assert (val != 0, ExcNotInitialized());
01329   Assert (factor !=0, ExcDivideByZero());
01330 
01331   const number factor_inv = 1. / factor;
01332 
01333   const unsigned int chunk_size = cols->get_chunk_size();
01334 
01335                                    // multiply all elements of the matrix with
01336                                    // the given factor. this includes the
01337                                    // padding elements in chunks that overlap
01338                                    // the boundaries of the actual matrix --
01339                                    // but since multiplication with a number
01340                                    // does not violate the invariant of
01341                                    // keeping these elements at zero nothing
01342                                    // can happen
01343   number             *val_ptr    = val;
01344   const number *const end_ptr    = val +
01345                                    cols->sparsity_pattern.n_nonzero_elements()
01346                                    *
01347                                    chunk_size * chunk_size;
01348 
01349   while (val_ptr != end_ptr)
01350     *val_ptr++ *= factor_inv;
01351 
01352   return *this;
01353 }
01354 
01355 
01356 
01357 template <typename number>
01358 inline
01359 number ChunkSparseMatrix<number>::operator () (const unsigned int i,
01360                                                const unsigned int j) const
01361 {
01362   Assert (cols != 0, ExcNotInitialized());
01363   AssertThrow (compute_location(i,j) != SparsityPattern::invalid_entry,
01364                ExcInvalidIndex(i,j));
01365   return val[compute_location(i,j)];
01366 }
01367 
01368 
01369 
01370 template <typename number>
01371 inline
01372 number ChunkSparseMatrix<number>::el (const unsigned int i,
01373                                       const unsigned int j) const
01374 {
01375   Assert (cols != 0, ExcNotInitialized());
01376   const unsigned int index = compute_location(i,j);
01377 
01378   if (index != ChunkSparsityPattern::invalid_entry)
01379     return val[index];
01380   else
01381     return 0;
01382 }
01383 
01384 
01385 
01386 template <typename number>
01387 inline
01388 number ChunkSparseMatrix<number>::diag_element (const unsigned int i) const
01389 {
01390   Assert (cols != 0, ExcNotInitialized());
01391   Assert (cols->optimize_diagonal(),  ExcNotQuadratic());
01392   Assert (i<m(), ExcInvalidIndex1(i));
01393 
01394                                    // Use that the first element in each row
01395                                    // of a quadratic matrix is the main
01396                                    // diagonal of the chunk sparsity pattern
01397   const unsigned int chunk_size = cols->get_chunk_size();
01398   return val[cols->sparsity_pattern.rowstart[i/chunk_size]
01399              *
01400              chunk_size * chunk_size
01401              +
01402              (i % chunk_size) * chunk_size
01403              +
01404              (i % chunk_size)];
01405 }
01406 
01407 
01408 
01409 template <typename number>
01410 inline
01411 number & ChunkSparseMatrix<number>::diag_element (const unsigned int i)
01412 {
01413   Assert (cols != 0, ExcNotInitialized());
01414   Assert (cols->optimize_diagonal(),  ExcNotQuadratic());
01415   Assert (i<m(), ExcInvalidIndex1(i));
01416 
01417                                    // Use that the first element in each row
01418                                    // of a quadratic matrix is the main
01419                                    // diagonal of the chunk sparsity pattern
01420   const unsigned int chunk_size = cols->get_chunk_size();
01421   return val[cols->sparsity_pattern.rowstart[i/chunk_size]
01422              *
01423              chunk_size * chunk_size
01424              +
01425              (i % chunk_size) * chunk_size
01426              +
01427              (i % chunk_size)];
01428 }
01429 
01430 
01431 
01432 template <typename number>
01433 template <typename ForwardIterator>
01434 void
01435 ChunkSparseMatrix<number>::copy_from (const ForwardIterator begin,
01436                                       const ForwardIterator end)
01437 {
01438   Assert (static_cast<unsigned int>(std::distance (begin, end)) == m(),
01439           ExcIteratorRange (std::distance (begin, end), m()));
01440 
01441                                    // for use in the inner loop, we
01442                                    // define a typedef to the type of
01443                                    // the inner iterators
01444   typedef typename std::iterator_traits<ForwardIterator>::value_type::const_iterator inner_iterator;
01445   unsigned int row=0;
01446   for (ForwardIterator i=begin; i!=end; ++i, ++row)
01447     {
01448       const inner_iterator end_of_row = i->end();
01449       for (inner_iterator j=i->begin(); j!=end_of_row; ++j)
01450                                          // write entries
01451         set (row, j->first, j->second);
01452     }
01453 }
01454 
01455 
01456 #endif // DOXYGEN
01457 
01458 
01459 /*----------------------------   chunk_sparse_matrix.h     ---------------------------*/
01460 
01461 DEAL_II_NAMESPACE_CLOSE
01462 
01463 #endif
01464 /*----------------------------   chunk_sparse_matrix.h     ---------------------------*/
01465 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Thu May 17 2012 20:04:40 by doxygen 1.7.3