00001
00002
00003
00004
00005
00006
00007
00008
00009
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
01191
01192 template <typename somenumber> friend class ChunkSparseMatrix;
01193 };
01194
01197 #ifndef DOXYGEN
01198
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
01255
01256
01257
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
01302
01303
01304
01305
01306
01307
01308
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
01336
01337
01338
01339
01340
01341
01342
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
01395
01396
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
01418
01419
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
01442
01443
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
01451 set (row, j->first, j->second);
01452 }
01453 }
01454
01455
01456 #endif // DOXYGEN
01457
01458
01459
01460
01461 DEAL_II_NAMESPACE_CLOSE
01462
01463 #endif
01464
01465