00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__sparsity_pattern_h
00013 #define __deal2__sparsity_pattern_h
00014
00015
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/base/exceptions.h>
00018 #include <deal.II/base/subscriptor.h>
00019 #include <boost/serialization/array.hpp>
00020 #include <boost/serialization/split_member.hpp>
00021
00022 #include <vector>
00023 #include <iostream>
00024
00025 DEAL_II_NAMESPACE_OPEN
00026
00027 class SparsityPattern;
00028 class ChunkSparsityPattern;
00029 template <typename number> class FullMatrix;
00030 template <typename number> class SparseMatrix;
00031 template <class VECTOR> class VectorSlice;
00032
00033 class CompressedSparsityPattern;
00034 class CompressedSetSparsityPattern;
00035 class CompressedSimpleSparsityPattern;
00036
00037
00038
00043 namespace internals
00044 {
00045 namespace SparsityPatternTools
00046 {
00055 unsigned int
00056 get_column_index_from_iterator (const unsigned int i);
00057
00067 template <typename value>
00068 unsigned int
00069 get_column_index_from_iterator (const std::pair<unsigned int, value> &i);
00070
00078 template <typename value>
00079 unsigned int
00080 get_column_index_from_iterator (const std::pair<const unsigned int, value> &i);
00081
00082 }
00083 }
00084
00085
00086
00087 namespace SparsityPatternIterators
00088 {
00089
00090 class Iterator;
00091
00105 class Accessor
00106 {
00107 public:
00111 Accessor (const SparsityPattern *matrix,
00112 const unsigned int row,
00113 const unsigned int index);
00114
00120 Accessor (const SparsityPattern *matrix);
00121
00129 unsigned int row () const;
00130
00138 unsigned int index () const;
00139
00147 unsigned int column () const;
00148
00166 inline bool is_valid_entry () const;
00167
00168
00175 bool operator == (const Accessor &) const;
00176
00177
00191 bool operator < (const Accessor &) const;
00192 protected:
00197 const SparsityPattern * sparsity_pattern;
00198
00202 unsigned int a_row;
00203
00207 unsigned int a_index;
00208
00213 void advance ();
00214
00218 friend class Iterator;
00219 };
00220
00221
00222
00227 class Iterator
00228 {
00229 public:
00235 Iterator (const SparsityPattern *sp,
00236 const unsigned int row,
00237 const unsigned int index);
00238
00242 Iterator& operator++ ();
00243
00247 Iterator operator++ (int);
00248
00252 const Accessor & operator* () const;
00253
00257 const Accessor * operator-> () const;
00258
00265 bool operator == (const Iterator&) const;
00266
00270 bool operator != (const Iterator&) const;
00271
00285 bool operator < (const Iterator&) const;
00286
00287 private:
00292 Accessor accessor;
00293 };
00294 }
00295
00296
00297
00307 class SparsityPattern : public Subscriptor
00308 {
00309 public:
00315 typedef
00316 SparsityPatternIterators::Iterator
00317 const_iterator;
00318
00324 typedef
00325 const unsigned int * row_iterator;
00326
00337 typedef
00338 SparsityPatternIterators::Iterator
00339 iterator;
00340
00341
00366 static const unsigned int invalid_entry = numbers::invalid_unsigned_int;
00367
00372
00383 SparsityPattern ();
00384
00416 SparsityPattern (const SparsityPattern &);
00417
00433 SparsityPattern (const unsigned int m,
00434 const unsigned int n,
00435 const unsigned int max_per_row,
00436 const bool optimize_diagonal = true);
00437
00456 SparsityPattern (const unsigned int m,
00457 const unsigned int n,
00458 const std::vector<unsigned int>& row_lengths,
00459 const bool optimize_diagonal = true);
00460
00474 SparsityPattern (const unsigned int n,
00475 const unsigned int max_per_row);
00476
00491 SparsityPattern (const unsigned int m,
00492 const std::vector<unsigned int>& row_lengths,
00493 const bool optimize_diagonal = true);
00494
00537 SparsityPattern (const SparsityPattern &original,
00538 const unsigned int max_per_row,
00539 const unsigned int extra_off_diagonals);
00540
00544 ~SparsityPattern ();
00545
00554 SparsityPattern & operator = (const SparsityPattern &);
00555
00567 void reinit (const unsigned int m,
00568 const unsigned int n,
00569 const unsigned int max_per_row,
00570 const bool optimize_diagonal = true);
00571
00597 void reinit (const unsigned int m,
00598 const unsigned int n,
00599 const std::vector<unsigned int> &row_lengths,
00600 const bool optimize_diagonal = true);
00601
00606 void reinit (const unsigned int m,
00607 const unsigned int n,
00608 const VectorSlice<const std::vector<unsigned int> > &row_lengths,
00609 const bool optimize_diagonal = true);
00610
00631 void compress ();
00632
00771 template <typename ForwardIterator>
00772 void copy_from (const unsigned int n_rows,
00773 const unsigned int n_cols,
00774 const ForwardIterator begin,
00775 const ForwardIterator end,
00776 const bool optimize_diagonal = true);
00777
00787 template <typename CompressedSparsityType>
00788 void copy_from (const CompressedSparsityType &csp,
00789 const bool optimize_diagonal = true);
00790
00802 template <typename number>
00803 void copy_from (const FullMatrix<number> &matrix,
00804 const bool optimize_diagonal = true);
00805
00817 void symmetrize ();
00818
00827 void add (const unsigned int i,
00828 const unsigned int j);
00829
00839 template <typename ForwardIterator>
00840 void add_entries (const unsigned int row,
00841 ForwardIterator begin,
00842 ForwardIterator end,
00843 const bool indices_are_sorted = false);
00844
00845
00849
00850
00857 inline iterator begin () const;
00858
00862 inline iterator end () const;
00863
00876 inline iterator begin (const unsigned int r) const;
00877
00889 inline iterator end (const unsigned int r) const;
00890
00903 inline row_iterator row_begin (const unsigned int r) const;
00904
00916 inline row_iterator row_end (const unsigned int r) const;
00917
00918
00922
00926 bool operator == (const SparsityPattern &) const;
00927
00934 bool empty () const;
00935
00944 unsigned int max_entries_per_row () const;
00945
00956 unsigned int bandwidth () const;
00957
00970 std::size_t n_nonzero_elements () const;
00971
00976 bool is_compressed () const;
00977
00983 inline unsigned int n_rows () const;
00984
00990 inline unsigned int n_cols () const;
00991
00995 unsigned int row_length (const unsigned int row) const;
00996
01031 bool optimize_diagonal () const;
01032
01051 bool stores_only_added_elements () const;
01052
01059 std::size_t memory_consumption () const;
01060
01061
01065
01094 unsigned int operator() (const unsigned int i,
01095 const unsigned int j) const;
01096
01117 std::pair<unsigned int, unsigned int>
01118 matrix_position (const unsigned int global_index) const;
01119
01124 bool exists (const unsigned int i,
01125 const unsigned int j) const;
01126
01137 unsigned int row_position(const unsigned int i,
01138 const unsigned int j) const;
01139
01159 unsigned int column_number (const unsigned int row,
01160 const unsigned int index) const;
01161
01162
01163
01167
01188 void block_write (std::ostream &out) const;
01189
01213 void block_read (std::istream &in);
01214
01224 void print (std::ostream &out) const;
01225
01250 void print_gnuplot (std::ostream &out) const;
01255 template <class Archive>
01256 void save (Archive & ar, const unsigned int version) const;
01257
01262 template <class Archive>
01263 void load (Archive & ar, const unsigned int version);
01264
01265
01269
01330 void partition (const unsigned int n_partitions,
01331 std::vector<unsigned int> &partition_indices) const;
01332
01333
01361 inline const std::size_t * get_rowstart_indices () const;
01362
01389 inline const unsigned int * get_column_numbers () const;
01390
01391 BOOST_SERIALIZATION_SPLIT_MEMBER()
01399 DeclException2 (ExcNotEnoughSpace,
01400 int, int,
01401 << "Upon entering a new entry to row " << arg1
01402 << ": there was no free entry any more. " << std::endl
01403 << "(Maximum number of entries for this row: "
01404 << arg2 << "; maybe the matrix is already compressed?)");
01411 DeclException0 (ExcNotCompressed);
01419 DeclException0 (ExcMatrixIsCompressed);
01423 DeclException0 (ExcInvalidConstructorCall);
01433 DeclException0 (ExcDiagonalNotOptimized);
01437 DeclException2 (ExcIteratorRange,
01438 int, int,
01439 << "The iterators denote a range of " << arg1
01440 << " elements, but the given number of rows was " << arg2);
01444 DeclException1 (ExcInvalidNumberOfPartitions,
01445 int,
01446 << "The number of partitions you gave is " << arg1
01447 << ", but must be greater than zero.");
01449 private:
01461 unsigned int max_dim;
01462
01467 unsigned int rows;
01468
01473 unsigned int cols;
01474
01482 std::size_t max_vec_len;
01483
01498 unsigned int max_row_length;
01499
01525 std::size_t *rowstart;
01526
01576 unsigned int *colnums;
01577
01583 bool compressed;
01584
01589 bool diagonal_optimized;
01590
01595 template <typename number> friend class SparseMatrix;
01596 template <typename number> friend class ChunkSparseMatrix;
01597
01598 friend class ChunkSparsityPattern;
01599 };
01600
01601
01603
01604
01605 #ifndef DOXYGEN
01606
01607
01608 namespace SparsityPatternIterators
01609 {
01610 inline
01611 Accessor::
01612 Accessor (const SparsityPattern *sparsity_pattern,
01613 const unsigned int r,
01614 const unsigned int i)
01615 :
01616 sparsity_pattern(sparsity_pattern),
01617 a_row(r),
01618 a_index(i)
01619 {}
01620
01621
01622 inline
01623 Accessor::
01624 Accessor (const SparsityPattern *sparsity_pattern)
01625 :
01626 sparsity_pattern(sparsity_pattern),
01627 a_row(sparsity_pattern->n_rows()),
01628 a_index(0)
01629 {}
01630
01631
01632 inline
01633 bool
01634 Accessor::is_valid_entry () const
01635 {
01636 return (sparsity_pattern
01637 ->get_column_numbers()[sparsity_pattern
01638 ->get_rowstart_indices()[a_row]+a_index]
01639 != SparsityPattern::invalid_entry);
01640 }
01641
01642
01643 inline
01644 unsigned int
01645 Accessor::row() const
01646 {
01647 Assert (is_valid_entry() == true, ExcInvalidIterator());
01648
01649 return a_row;
01650 }
01651
01652
01653 inline
01654 unsigned int
01655 Accessor::column() const
01656 {
01657 Assert (is_valid_entry() == true, ExcInvalidIterator());
01658
01659 return (sparsity_pattern
01660 ->get_column_numbers()[sparsity_pattern
01661 ->get_rowstart_indices()[a_row]+a_index]);
01662 }
01663
01664
01665 inline
01666 unsigned int
01667 Accessor::index() const
01668 {
01669 Assert (is_valid_entry() == true, ExcInvalidIterator());
01670
01671 return a_index;
01672 }
01673
01674
01675
01676
01677 inline
01678 bool
01679 Accessor::operator == (const Accessor& other) const
01680 {
01681 return (sparsity_pattern == other.sparsity_pattern &&
01682 a_row == other.a_row &&
01683 a_index == other.a_index);
01684 }
01685
01686
01687
01688 inline
01689 bool
01690 Accessor::operator < (const Accessor& other) const
01691 {
01692 Assert (sparsity_pattern == other.sparsity_pattern,
01693 ExcInternalError());
01694
01695 return (a_row < other.a_row ||
01696 (a_row == other.a_row &&
01697 a_index < other.a_index));
01698 }
01699
01700
01701 inline
01702 void
01703 Accessor::advance ()
01704 {
01705 Assert (a_row < sparsity_pattern->n_rows(), ExcIteratorPastEnd());
01706
01707 ++a_index;
01708
01709
01710
01711
01712 while (a_index >= sparsity_pattern->row_length(a_row))
01713 {
01714 a_index = 0;
01715 ++a_row;
01716
01717
01718
01719 if (a_row == sparsity_pattern->n_rows())
01720 break;
01721 }
01722 }
01723
01724
01725
01726 inline
01727 Iterator::Iterator (const SparsityPattern *sparsity_pattern,
01728 const unsigned int r,
01729 const unsigned int i)
01730 :
01731 accessor(sparsity_pattern, r, i)
01732 {}
01733
01734
01735
01736 inline
01737 Iterator &
01738 Iterator::operator++ ()
01739 {
01740 accessor.advance ();
01741 return *this;
01742 }
01743
01744
01745
01746 inline
01747 Iterator
01748 Iterator::operator++ (int)
01749 {
01750 const Iterator iter = *this;
01751 accessor.advance ();
01752 return iter;
01753 }
01754
01755
01756
01757 inline
01758 const Accessor &
01759 Iterator::operator* () const
01760 {
01761 return accessor;
01762 }
01763
01764
01765
01766 inline
01767 const Accessor *
01768 Iterator::operator-> () const
01769 {
01770 return &accessor;
01771 }
01772
01773
01774 inline
01775 bool
01776 Iterator::operator == (const Iterator& other) const
01777 {
01778 return (accessor == other.accessor);
01779 }
01780
01781
01782
01783 inline
01784 bool
01785 Iterator::operator != (const Iterator& other) const
01786 {
01787 return ! (*this == other);
01788 }
01789
01790
01791 inline
01792 bool
01793 Iterator::operator < (const Iterator& other) const
01794 {
01795 return accessor < other.accessor;
01796 }
01797
01798 }
01799
01800
01801
01802 inline
01803 SparsityPattern::iterator
01804 SparsityPattern::begin () const
01805 {
01806
01807
01808 for (unsigned int r=0; r<n_rows(); ++r)
01809 if (row_length(r) > 0)
01810 return iterator(this, r, 0);
01811
01812
01813
01814
01815 return end();
01816 }
01817
01818
01819 inline
01820 SparsityPattern::iterator
01821 SparsityPattern::end () const
01822 {
01823 return iterator(this, n_rows(), 0);
01824 }
01825
01826
01827
01828 inline
01829 SparsityPattern::iterator
01830 SparsityPattern::begin (const unsigned int r) const
01831 {
01832 Assert (r<n_rows(), ExcIndexRange(r,0,n_rows()));
01833
01834 if (row_length(r) > 0)
01835 return iterator(this, r, 0);
01836 else
01837 return end (r);
01838 }
01839
01840
01841
01842 inline
01843 SparsityPattern::iterator
01844 SparsityPattern::end (const unsigned int r) const
01845 {
01846 Assert (r<n_rows(), ExcIndexRange(r,0,n_rows()));
01847
01848
01849
01850
01851 for (unsigned int i=r+1; i<n_rows(); ++i)
01852 if (row_length(i) > 0)
01853 return iterator(this, i, 0);
01854
01855
01856
01857 return end();
01858 }
01859
01860
01861
01862 inline
01863 SparsityPattern::row_iterator
01864 SparsityPattern::row_begin (const unsigned int r) const
01865 {
01866 Assert (r<n_rows(), ExcIndexRange(r,0,n_rows()));
01867 return &colnums[rowstart[r]];
01868 }
01869
01870
01871
01872 inline
01873 SparsityPattern::row_iterator
01874 SparsityPattern::row_end (const unsigned int r) const
01875 {
01876 Assert (r<n_rows(), ExcIndexRange(r,0,n_rows()));
01877 return &colnums[rowstart[r+1]];
01878 }
01879
01880
01881
01882 inline
01883 unsigned int
01884 SparsityPattern::n_rows () const
01885 {
01886 return rows;
01887 }
01888
01889
01890 inline
01891 unsigned int
01892 SparsityPattern::n_cols () const
01893 {
01894 return cols;
01895 }
01896
01897
01898 inline
01899 bool
01900 SparsityPattern::is_compressed () const
01901 {
01902 return compressed;
01903 }
01904
01905
01906 inline
01907 bool
01908 SparsityPattern::optimize_diagonal () const
01909 {
01910 return diagonal_optimized;
01911 }
01912
01913
01914 inline
01915 bool
01916 SparsityPattern::stores_only_added_elements () const
01917 {
01918 if ((diagonal_optimized == true)
01919 &&
01920 (n_cols() == n_rows()))
01921 return false;
01922 else
01923 return true;
01924 }
01925
01926
01927 inline
01928 const std::size_t *
01929 SparsityPattern::get_rowstart_indices () const
01930 {
01931 return rowstart;
01932 }
01933
01934
01935 inline
01936 const unsigned int *
01937 SparsityPattern::get_column_numbers () const
01938 {
01939 return colnums;
01940 }
01941
01942
01943
01944 inline
01945 unsigned int
01946 SparsityPattern::row_length (const unsigned int row) const
01947 {
01948 Assert(row<rows, ExcIndexRange(row,0,rows));
01949 return rowstart[row+1]-rowstart[row];
01950 }
01951
01952
01953
01954 inline
01955 unsigned int
01956 SparsityPattern::column_number (const unsigned int row,
01957 const unsigned int index) const
01958 {
01959 Assert(row<rows, ExcIndexRange(row,0,rows));
01960 Assert(index<row_length(row), ExcIndexRange(index,0,row_length(row)));
01961
01962 return colnums[rowstart[row]+index];
01963 }
01964
01965
01966 inline
01967 std::size_t
01968 SparsityPattern::n_nonzero_elements () const
01969 {
01970 Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
01971 Assert (compressed, ExcNotCompressed());
01972 return rowstart[rows]-rowstart[0];
01973 }
01974
01975
01976
01977 template <class Archive>
01978 inline
01979 void
01980 SparsityPattern::save (Archive & ar, const unsigned int) const
01981 {
01982
01983
01984 ar & static_cast<const Subscriptor &>(*this);
01985
01986 ar & max_dim & rows & cols & max_vec_len & max_row_length & compressed & diagonal_optimized;
01987
01988 ar & boost::serialization::make_array(rowstart, max_dim + 1);
01989 ar & boost::serialization::make_array(colnums, max_vec_len);
01990 }
01991
01992
01993
01994 template <class Archive>
01995 inline
01996 void
01997 SparsityPattern::load (Archive & ar, const unsigned int)
01998 {
01999
02000
02001 ar & static_cast<Subscriptor &>(*this);
02002
02003 ar & max_dim & rows & cols & max_vec_len & max_row_length & compressed & diagonal_optimized;
02004
02005 rowstart = new std::size_t [max_dim + 1];
02006 colnums = new unsigned int [max_vec_len];
02007
02008 ar & boost::serialization::make_array(rowstart, max_dim + 1);
02009 ar & boost::serialization::make_array(colnums, max_vec_len);
02010 }
02011
02012
02013
02014 inline
02015 bool
02016 SparsityPattern::operator == (const SparsityPattern& sp2) const
02017 {
02018
02019
02020
02021
02022
02023
02024
02025
02026
02027 if (rows != sp2.rows ||
02028 cols != sp2.cols ||
02029 compressed != sp2.compressed ||
02030 diagonal_optimized != sp2.diagonal_optimized)
02031 return false;
02032
02033 for (unsigned int i = 0; i < rows+1; ++i)
02034 if (rowstart[i] != sp2.rowstart[i])
02035 return false;
02036
02037 for (unsigned int i = 0; i < rowstart[rows]; ++i)
02038 if (colnums[i] != sp2.colnums[i])
02039 return false;
02040
02041 return true;
02042 }
02043
02044
02045
02046 namespace internal
02047 {
02048 namespace SparsityPatternTools
02049 {
02050 inline
02051 unsigned int
02052 get_column_index_from_iterator (const unsigned int i)
02053 {
02054 return i;
02055 }
02056
02057
02058
02059 template <typename value>
02060 inline
02061 unsigned int
02062 get_column_index_from_iterator (const std::pair<unsigned int, value> &i)
02063 {
02064 return i.first;
02065 }
02066
02067
02068
02069 template <typename value>
02070 inline
02071 unsigned int
02072 get_column_index_from_iterator (const std::pair<const unsigned int, value> &i)
02073 {
02074 return i.first;
02075 }
02076 }
02077 }
02078
02079
02080
02081 template <typename ForwardIterator>
02082 void
02083 SparsityPattern::copy_from (const unsigned int n_rows,
02084 const unsigned int n_cols,
02085 const ForwardIterator begin,
02086 const ForwardIterator end,
02087 const bool optimize_diag)
02088 {
02089 Assert (static_cast<unsigned int>(std::distance (begin, end)) == n_rows,
02090 ExcIteratorRange (std::distance (begin, end), n_rows));
02091
02092
02093
02094
02095
02096
02097
02098
02099
02100
02101
02102 const bool is_square = optimize_diag && (n_rows == n_cols);
02103 std::vector<unsigned int> row_lengths;
02104 row_lengths.reserve(n_rows);
02105 for (ForwardIterator i=begin; i!=end; ++i)
02106 row_lengths.push_back (std::distance (i->begin(), i->end())
02107 +
02108 (is_square ? 1 : 0));
02109 reinit (n_rows, n_cols, row_lengths, is_square);
02110
02111
02112
02113
02114
02115
02116
02117
02118
02119
02120 unsigned int row = 0;
02121 typedef typename std::iterator_traits<ForwardIterator>::value_type::const_iterator inner_iterator;
02122 for (ForwardIterator i=begin; i!=end; ++i, ++row)
02123 {
02124 unsigned int *cols = &colnums[rowstart[row]] + (is_square ? 1 : 0);
02125 const inner_iterator end_of_row = i->end();
02126 for (inner_iterator j=i->begin(); j!=end_of_row; ++j)
02127 {
02128 const unsigned int col
02129 = internal::SparsityPatternTools::get_column_index_from_iterator(*j);
02130 Assert (col < n_cols, ExcIndexRange(col,0,n_cols));
02131
02132 if ((col!=row) || !is_square)
02133 *cols++ = col;
02134 };
02135 };
02136
02137
02138
02139
02140 compress ();
02141 }
02142
02143
02144 #endif // DOXYGEN
02145
02146 DEAL_II_NAMESPACE_CLOSE
02147
02148 #endif
02149