00001
00002
00003
00004
00005
00006
00007
00008
00009
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
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
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
02126
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
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
02165
02166
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
02179
02180
02181
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
02268
02269
02270
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
02421
02422
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
02437
02438
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
02495
02496
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
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
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
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
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
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
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
02824
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
02830
02831
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
02851
02852 for (unsigned int r=0; r<m(); ++r)
02853 if (cols->row_length(r) > 0)
02854 return iterator(this, r, 0);
02855
02856
02857
02858
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
02895
02896
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
02902
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
02931
02932
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
02938
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
02973 }
02974
02975
02976
02977 template <typename number>
02978 inline
02979 void
02980 SparseMatrix<number>::
02981 prepare_set()
02982 {
02983
02984 }
02985
02986 #endif // DOXYGEN
02987
02988
02989
02990
02991 DEAL_II_NAMESPACE_CLOSE
02992
02993 #endif
02994
02995