00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__full_matrix_h
00013 #define __deal2__full_matrix_h
00014
00015
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/base/numbers.h>
00018 #include <deal.II/base/table.h>
00019 #include <deal.II/lac/exceptions.h>
00020 #include <deal.II/lac/identity_matrix.h>
00021 #include <deal.II/base/tensor.h>
00022
00023 #include <vector>
00024 #include <iomanip>
00025 #include <cstring>
00026
00027 DEAL_II_NAMESPACE_OPEN
00028
00029
00030
00031 template <typename number> class Vector;
00032 template <typename number> class LAPACKFullMatrix;
00033
00034
00059 template<typename number>
00060 class FullMatrix : public Table<2,number>
00061 {
00062 public:
00067 typedef number value_type;
00068
00069
00089 typedef typename numbers::NumberTraits<number>::real_type real_type;
00090
00091
00092 class const_iterator;
00093
00097 class Accessor
00098 {
00099 public:
00106 Accessor (const FullMatrix<number> *matrix,
00107 const unsigned int row,
00108 const unsigned int col);
00109
00115 unsigned int row() const;
00116
00122 unsigned int column() const;
00123
00127 number value() const;
00128
00129 protected:
00133 const FullMatrix<number>* matrix;
00134
00138 unsigned int a_row;
00139
00143 unsigned short a_col;
00144
00145
00146
00147
00148
00149 friend class const_iterator;
00150 };
00151
00155 class const_iterator
00156 {
00157 public:
00161 const_iterator(const FullMatrix<number> *matrix,
00162 const unsigned int row,
00163 const unsigned int col);
00164
00168 const_iterator& operator++ ();
00169
00173 const_iterator& operator++ (int);
00174
00178 const Accessor& operator* () const;
00179
00183 const Accessor* operator-> () const;
00184
00191 bool operator == (const const_iterator&) const;
00195 bool operator != (const const_iterator&) const;
00196
00204 bool operator < (const const_iterator&) const;
00205
00211 bool operator > (const const_iterator&) const;
00212
00213 private:
00218 Accessor accessor;
00219 };
00225
00240 explicit FullMatrix (const unsigned int n = 0);
00241
00247 FullMatrix (const unsigned int rows,
00248 const unsigned int cols);
00249
00269 FullMatrix (const FullMatrix&);
00270
00277 FullMatrix (const unsigned int rows,
00278 const unsigned int cols,
00279 const number* entries);
00280
00293 explicit FullMatrix (const IdentityMatrix &id);
00308 FullMatrix<number> &
00309 operator = (const FullMatrix<number>&);
00310
00314 template<typename number2>
00315 FullMatrix<number> &
00316 operator = (const FullMatrix<number2>&);
00317
00328 FullMatrix<number> &
00329 operator = (const number d);
00330
00342 FullMatrix<number> &
00343 operator = (const IdentityMatrix &id);
00344
00351 template <typename number2>
00352 FullMatrix<number> &
00353 operator = (const LAPACKFullMatrix<number2>&);
00354
00355
00364 template <class MATRIX>
00365 void copy_from (const MATRIX&);
00366
00375 template <class MATRIX>
00376 void copy_transposed (const MATRIX&);
00377
00391 template <int dim>
00392 void
00393 copy_from (const Tensor<2,dim> &T,
00394 const unsigned int src_r_i=0,
00395 const unsigned int src_r_j=dim-1,
00396 const unsigned int src_c_i=0,
00397 const unsigned int src_c_j=dim-1,
00398 const unsigned int dst_r=0,
00399 const unsigned int dst_c=0);
00400
00414 template <int dim>
00415 void
00416 copy_to(Tensor<2,dim> &T,
00417 const unsigned int src_r_i=0,
00418 const unsigned int src_r_j=dim-1,
00419 const unsigned int src_c_i=0,
00420 const unsigned int src_c_j=dim-1,
00421 const unsigned int dst_r=0,
00422 const unsigned int dst_c=0) const;
00423
00436 template <typename MatrixType>
00437 void extract_submatrix_from (const MatrixType &matrix,
00438 const std::vector<unsigned int> &row_index_set,
00439 const std::vector<unsigned int> &column_index_set);
00440
00453 template <typename MatrixType>
00454 void
00455 scatter_matrix_to (const std::vector<unsigned int> &row_index_set,
00456 const std::vector<unsigned int> &column_index_set,
00457 MatrixType &matrix) const;
00458
00477 template<typename number2>
00478 void fill (const FullMatrix<number2> &src,
00479 const unsigned int dst_offset_i = 0,
00480 const unsigned int dst_offset_j = 0,
00481 const unsigned int src_offset_i = 0,
00482 const unsigned int src_offset_j = 0);
00483
00484
00489 template<typename number2>
00490 void fill (const number2*);
00491
00510 template<typename number2>
00511 void fill_permutation (const FullMatrix<number2> &src,
00512 const std::vector<unsigned int> &p_rows,
00513 const std::vector<unsigned int> &p_cols);
00514
00525 void set (const unsigned int i,
00526 const unsigned int j,
00527 const number value);
00547 bool operator == (const FullMatrix<number> &) const;
00548
00554 unsigned int m () const;
00555
00561 unsigned int n () const;
00562
00573 bool all_zero () const;
00574
00602 template<typename number2>
00603 number2 matrix_norm_square (const Vector<number2> &v) const;
00604
00620 template<typename number2>
00621 number2 matrix_scalar_product (const Vector<number2> &u,
00622 const Vector<number2> &v) const;
00623
00632 real_type l1_norm () const;
00633
00641 real_type linfty_norm () const;
00642
00655 real_type frobenius_norm () const;
00656
00670 real_type relative_symmetry_norm2 () const;
00671
00682 number determinant () const;
00683
00692 number trace () const;
00693
00698 template <class STREAM>
00699 void print (STREAM &s,
00700 const unsigned int width=5,
00701 const unsigned int precision=2) const;
00702
00744 void print_formatted (std::ostream &out,
00745 const unsigned int precision=3,
00746 const bool scientific = true,
00747 const unsigned int width = 0,
00748 const char *zero_string = " ",
00749 const double denominator = 1.,
00750 const double threshold = 0.) const;
00751
00757 std::size_t memory_consumption () const;
00758
00760
00761
00762
00767 const_iterator begin () const;
00768
00772 const_iterator end () const;
00773
00778 const_iterator begin (const unsigned int r) const;
00779
00783 const_iterator end (const unsigned int r) const;
00784
00786
00787
00788
00793 FullMatrix & operator *= (const number factor);
00794
00799 FullMatrix & operator /= (const number factor);
00800
00813 template<typename number2>
00814 void add (const number a,
00815 const FullMatrix<number2> &A);
00816
00830 template<typename number2>
00831 void add (const number a,
00832 const FullMatrix<number2> &A,
00833 const number b,
00834 const FullMatrix<number2> &B);
00835
00849 template<typename number2>
00850 void add (const number a,
00851 const FullMatrix<number2> &A,
00852 const number b,
00853 const FullMatrix<number2> &B,
00854 const number c,
00855 const FullMatrix<number2> &C);
00856
00874 template<typename number2>
00875 void add (const FullMatrix<number2> &src,
00876 const number factor,
00877 const unsigned int dst_offset_i = 0,
00878 const unsigned int dst_offset_j = 0,
00879 const unsigned int src_offset_i = 0,
00880 const unsigned int src_offset_j = 0);
00881
00889 template<typename number2>
00890 void Tadd (const number s,
00891 const FullMatrix<number2> &B);
00892
00915 template<typename number2>
00916 void Tadd (const FullMatrix<number2> &src,
00917 const number factor,
00918 const unsigned int dst_offset_i = 0,
00919 const unsigned int dst_offset_j = 0,
00920 const unsigned int src_offset_i = 0,
00921 const unsigned int src_offset_j = 0);
00922
00927 void add(const unsigned int row, const unsigned int column, const number value);
00928
00945 template <typename number2>
00946 void add (const unsigned int row,
00947 const unsigned int n_cols,
00948 const unsigned int *col_indices,
00949 const number2 *values,
00950 const bool elide_zero_values = true,
00951 const bool col_indices_are_sorted = false);
00952
00958 void add_row (const unsigned int i,
00959 const number s,
00960 const unsigned int j);
00961
00967 void add_row (const unsigned int i,
00968 const number s, const unsigned int j,
00969 const number t, const unsigned int k);
00970
00975 void add_col (const unsigned int i,
00976 const number s,
00977 const unsigned int j);
00978
00984 void add_col (const unsigned int i,
00985 const number s, const unsigned int j,
00986 const number t, const unsigned int k);
00987
00993 void swap_row (const unsigned int i,
00994 const unsigned int j);
00995
01001 void swap_col (const unsigned int i,
01002 const unsigned int j);
01003
01010 void diagadd (const number s);
01011
01016 template<typename number2>
01017 void equ (const number a,
01018 const FullMatrix<number2> &A);
01019
01024 template<typename number2>
01025 void equ (const number a,
01026 const FullMatrix<number2> &A,
01027 const number b,
01028 const FullMatrix<number2> &B);
01029
01034 template<typename number2>
01035 void equ (const number a,
01036 const FullMatrix<number2> &A,
01037 const number b,
01038 const FullMatrix<number2> &B,
01039 const number c,
01040 const FullMatrix<number2> &C);
01041
01052 void symmetrize ();
01053
01077 void gauss_jordan ();
01078
01089 template <typename number2>
01090 void invert (const FullMatrix<number2> &M);
01091
01102 template <typename number2>
01103 void cholesky (const FullMatrix<number2> &A);
01104
01110 template <typename number2>
01111 void outer_product (const Vector<number2> &V,
01112 const Vector<number2> &W);
01113
01120 template <typename number2>
01121 void left_invert (const FullMatrix<number2> &M);
01122
01129 template <typename number2>
01130 void right_invert (const FullMatrix<number2> &M);
01131
01133
01134
01135
01162 template<typename number2>
01163 void mmult (FullMatrix<number2> &C,
01164 const FullMatrix<number2> &B,
01165 const bool adding=false) const;
01166
01194 template<typename number2>
01195 void Tmmult (FullMatrix<number2> &C,
01196 const FullMatrix<number2> &B,
01197 const bool adding=false) const;
01198
01226 template<typename number2>
01227 void mTmult (FullMatrix<number2> &C,
01228 const FullMatrix<number2> &B,
01229 const bool adding=false) const;
01230
01259 template<typename number2>
01260 void TmTmult (FullMatrix<number2> &C,
01261 const FullMatrix<number2> &B,
01262 const bool adding=false) const;
01263
01284 void triple_product(const FullMatrix<number>& A,
01285 const FullMatrix<number>& B,
01286 const FullMatrix<number>& D,
01287 const bool transpose_B = false,
01288 const bool transpose_D = false,
01289 const number scaling = number(1.));
01290
01308 template<typename number2>
01309 void vmult (Vector<number2> &w,
01310 const Vector<number2> &v,
01311 const bool adding=false) const;
01312
01320 template<typename number2>
01321 void vmult_add (Vector<number2> &w,
01322 const Vector<number2> &v) const;
01323
01343 template<typename number2>
01344 void Tvmult (Vector<number2> &w,
01345 const Vector<number2> &v,
01346 const bool adding=false) const;
01347
01356 template<typename number2>
01357 void Tvmult_add (Vector<number2> &w,
01358 const Vector<number2> &v) const;
01359
01370 template <typename somenumber>
01371 void precondition_Jacobi (Vector<somenumber> &dst,
01372 const Vector<somenumber> &src,
01373 const number omega = 1.) const;
01374
01384 template<typename number2, typename number3>
01385 number residual (Vector<number2> &dst,
01386 const Vector<number2> &x,
01387 const Vector<number3> &b) const;
01388
01406 template<typename number2>
01407 void forward (Vector<number2> &dst,
01408 const Vector<number2> &src) const;
01409
01420 template<typename number2>
01421 void backward (Vector<number2> &dst,
01422 const Vector<number2> &src) const;
01423
01425
01432 DeclException0 (ExcEmptyMatrix);
01433
01437 DeclException1 (ExcNotRegular,
01438 number,
01439 << "The maximal pivot is " << arg1
01440 << ", which is below the threshold. The matrix may be singular.");
01444 DeclException3 (ExcInvalidDestination,
01445 int, int, int,
01446 << "Target region not in matrix: size in this direction="
01447 << arg1 << ", size of new matrix=" << arg2
01448 << ", offset=" << arg3);
01452 DeclException0 (ExcSourceEqualsDestination);
01456 DeclException0 (ExcMatrixNotPositiveDefinite);
01458
01459 friend class Accessor;
01460 };
01461
01464 #ifndef DOXYGEN
01465
01466
01467
01468
01469
01470 template <typename number>
01471 inline
01472 unsigned int
01473 FullMatrix<number>::m() const
01474 {
01475 return this->n_rows();
01476 }
01477
01478
01479
01480 template <typename number>
01481 inline
01482 unsigned int
01483 FullMatrix<number>::n() const
01484 {
01485 return this->n_cols();
01486 }
01487
01488
01489
01490 template <typename number>
01491 FullMatrix<number> &
01492 FullMatrix<number>::operator = (const number d)
01493 {
01494 Assert (d==number(0), ExcScalarAssignmentOnlyForZeroValue());
01495
01496 if (this->n_elements() != 0)
01497 memset (&this->values[0], 0, this->n_elements()*sizeof(number));
01498
01499 return *this;
01500 }
01501
01502
01503
01504 template <typename number>
01505 template <typename number2>
01506 inline
01507 void FullMatrix<number>::fill (const number2* src)
01508 {
01509 Table<2,number>::fill(src);
01510 }
01511
01512
01513
01514 template <typename number>
01515 template <class MATRIX>
01516 void
01517 FullMatrix<number>::copy_from (const MATRIX& M)
01518 {
01519 this->reinit (M.m(), M.n());
01520 const typename MATRIX::const_iterator end = M.end();
01521 for (typename MATRIX::const_iterator entry = M.begin();
01522 entry != end; ++entry)
01523 this->el(entry->row(), entry->column()) = entry->value();
01524 }
01525
01526
01527
01528 template <typename number>
01529 template <class MATRIX>
01530 void
01531 FullMatrix<number>::copy_transposed (const MATRIX& M)
01532 {
01533 this->reinit (M.n(), M.m());
01534 const typename MATRIX::const_iterator end = M.end();
01535 for (typename MATRIX::const_iterator entry = M.begin();
01536 entry != end; ++entry)
01537 this->el(entry->column(), entry->row()) = entry->value();
01538 }
01539
01540
01541
01542 template <typename number>
01543 template <typename MatrixType>
01544 inline
01545 void
01546 FullMatrix<number>::extract_submatrix_from (const MatrixType &matrix,
01547 const std::vector<unsigned int> &row_index_set,
01548 const std::vector<unsigned int> &column_index_set)
01549 {
01550 AssertDimension(row_index_set.size(), this->n_rows());
01551 AssertDimension(column_index_set.size(), this->n_cols());
01552
01553 const unsigned int n_rows_submatrix = row_index_set.size();
01554 const unsigned int n_cols_submatrix = column_index_set.size();
01555
01556 for (unsigned int sub_row = 0; sub_row < n_rows_submatrix; ++sub_row)
01557 for (unsigned int sub_col = 0; sub_col < n_cols_submatrix; ++sub_col)
01558 (*this)(sub_row, sub_col) = matrix.el(row_index_set[sub_row], column_index_set[sub_col]);
01559 }
01560
01561
01562
01563 template <typename number>
01564 template <typename MatrixType>
01565 inline
01566 void
01567 FullMatrix<number>::scatter_matrix_to (const std::vector<unsigned int> &row_index_set,
01568 const std::vector<unsigned int> &column_index_set,
01569 MatrixType &matrix) const
01570 {
01571 AssertDimension(row_index_set.size(), this->n_rows());
01572 AssertDimension(column_index_set.size(), this->n_cols());
01573
01574 const unsigned int n_rows_submatrix = row_index_set.size();
01575 const unsigned int n_cols_submatrix = column_index_set.size();
01576
01577 for (unsigned int sub_row = 0; sub_row < n_rows_submatrix; ++sub_row)
01578 for (unsigned int sub_col = 0; sub_col < n_cols_submatrix; ++sub_col)
01579 matrix.set(row_index_set[sub_row],
01580 column_index_set[sub_col],
01581 (*this)(sub_row, sub_col));
01582 }
01583
01584
01585 template <typename number>
01586 inline
01587 void
01588 FullMatrix<number>::set (const unsigned int i,
01589 const unsigned int j,
01590 const number value)
01591 {
01592 (*this)(i,j) = value;
01593 }
01594
01595
01596
01597 template <typename number>
01598 template<typename number2>
01599 void
01600 FullMatrix<number>::vmult_add (Vector<number2> &w,
01601 const Vector<number2> &v) const
01602 {
01603 vmult(w, v, true);
01604 }
01605
01606
01607 template <typename number>
01608 template<typename number2>
01609 void
01610 FullMatrix<number>::Tvmult_add (Vector<number2> &w,
01611 const Vector<number2> &v) const
01612 {
01613 Tvmult(w, v, true);
01614 }
01615
01616
01617
01618
01619
01620 template <typename number>
01621 inline
01622 FullMatrix<number>::Accessor::
01623 Accessor (const FullMatrix<number>* matrix,
01624 const unsigned int r,
01625 const unsigned int c)
01626 :
01627 matrix(matrix),
01628 a_row(r),
01629 a_col(c)
01630 {}
01631
01632
01633 template <typename number>
01634 inline
01635 unsigned int
01636 FullMatrix<number>::Accessor::row() const
01637 {
01638 return a_row;
01639 }
01640
01641
01642 template <typename number>
01643 inline
01644 unsigned int
01645 FullMatrix<number>::Accessor::column() const
01646 {
01647 return a_col;
01648 }
01649
01650
01651 template <typename number>
01652 inline
01653 number
01654 FullMatrix<number>::Accessor::value() const
01655 {
01656 Assert (numbers::is_finite( matrix->el(a_row, a_col) ), ExcNumberNotFinite());
01657 return matrix->el(a_row, a_col);
01658 }
01659
01660
01661 template <typename number>
01662 inline
01663 FullMatrix<number>::const_iterator::
01664 const_iterator(const FullMatrix<number> *matrix,
01665 const unsigned int r,
01666 const unsigned int c)
01667 :
01668 accessor(matrix, r, c)
01669 {}
01670
01671
01672 template <typename number>
01673 inline
01674 typename FullMatrix<number>::const_iterator &
01675 FullMatrix<number>::const_iterator::operator++ ()
01676 {
01677 Assert (accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
01678
01679 ++accessor.a_col;
01680 if (accessor.a_col >= accessor.matrix->n())
01681 {
01682 accessor.a_col = 0;
01683 accessor.a_row++;
01684 }
01685 return *this;
01686 }
01687
01688
01689 template <typename number>
01690 inline
01691 const typename FullMatrix<number>::Accessor &
01692 FullMatrix<number>::const_iterator::operator* () const
01693 {
01694 return accessor;
01695 }
01696
01697
01698 template <typename number>
01699 inline
01700 const typename FullMatrix<number>::Accessor *
01701 FullMatrix<number>::const_iterator::operator-> () const
01702 {
01703 return &accessor;
01704 }
01705
01706
01707 template <typename number>
01708 inline
01709 bool
01710 FullMatrix<number>::const_iterator::
01711 operator == (const const_iterator& other) const
01712 {
01713 return (accessor.row() == other.accessor.row() &&
01714 accessor.column() == other.accessor.column());
01715 }
01716
01717
01718 template <typename number>
01719 inline
01720 bool
01721 FullMatrix<number>::const_iterator::
01722 operator != (const const_iterator& other) const
01723 {
01724 return ! (*this == other);
01725 }
01726
01727
01728 template <typename number>
01729 inline
01730 bool
01731 FullMatrix<number>::const_iterator::
01732 operator < (const const_iterator& other) const
01733 {
01734 return (accessor.row() < other.accessor.row() ||
01735 (accessor.row() == other.accessor.row() &&
01736 accessor.column() < other.accessor.column()));
01737 }
01738
01739
01740 template <typename number>
01741 inline
01742 bool
01743 FullMatrix<number>::const_iterator::
01744 operator > (const const_iterator& other) const
01745 {
01746 return (other < *this);
01747 }
01748
01749
01750 template <typename number>
01751 inline
01752 typename FullMatrix<number>::const_iterator
01753 FullMatrix<number>::begin () const
01754 {
01755 return const_iterator(this, 0, 0);
01756 }
01757
01758
01759 template <typename number>
01760 inline
01761 typename FullMatrix<number>::const_iterator
01762 FullMatrix<number>::end () const
01763 {
01764 return const_iterator(this, m(), 0);
01765 }
01766
01767
01768 template <typename number>
01769 inline
01770 typename FullMatrix<number>::const_iterator
01771 FullMatrix<number>::begin (const unsigned int r) const
01772 {
01773 Assert (r<m(), ExcIndexRange(r,0,m()));
01774 return const_iterator(this, r, 0);
01775 }
01776
01777
01778
01779 template <typename number>
01780 inline
01781 typename FullMatrix<number>::const_iterator
01782 FullMatrix<number>::end (const unsigned int r) const
01783 {
01784 Assert (r<m(), ExcIndexRange(r,0,m()));
01785 return const_iterator(this, r+1, 0);
01786 }
01787
01788
01789
01790 template <typename number>
01791 inline
01792 void
01793 FullMatrix<number>::add (const unsigned int r, const unsigned int c, const number v)
01794 {
01795 AssertIndexRange(r, this->m());
01796 AssertIndexRange(c, this->n());
01797
01798 this->operator()(r,c) += v;
01799 }
01800
01801
01802
01803 template <typename number>
01804 template <typename number2>
01805 inline
01806 void
01807 FullMatrix<number>::add (const unsigned int row,
01808 const unsigned int n_cols,
01809 const unsigned int *col_indices,
01810 const number2 *values,
01811 const bool,
01812 const bool)
01813 {
01814 AssertIndexRange(row, this->m());
01815 for (unsigned int col=0; col<n_cols; ++col)
01816 {
01817 AssertIndexRange(col_indices[col], this->n());
01818 this->operator()(row,col_indices[col]) += values[col];
01819 }
01820 }
01821
01822
01823 template <typename number>
01824 template <class STREAM>
01825 inline
01826 void
01827 FullMatrix<number>::print (STREAM &s,
01828 const unsigned int w,
01829 const unsigned int p) const
01830 {
01831 Assert (!this->empty(), ExcEmptyMatrix());
01832
01833 for (unsigned int i=0; i<this->m(); ++i)
01834 {
01835 for (unsigned int j=0; j<this->n(); ++j)
01836 s << std::setw(w) << std::setprecision(p) << this->el(i,j);
01837 s << std::endl;
01838 }
01839 }
01840
01841
01842 #endif // DOXYGEN
01843
01844 DEAL_II_NAMESPACE_CLOSE
01845
01846 #endif
01847
01848