Reference documentation for deal.II version Git 5546c7f774 2020-10-25 14:46:39 +0100
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
lapack_full_matrix.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2005 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_lapack_full_matrix_h
17 #define dealii_lapack_full_matrix_h
18 
19 
20 #include <deal.II/base/config.h>
21 
23 #include <deal.II/base/table.h>
25 
28 
29 #include <complex>
30 #include <memory>
31 #include <vector>
32 
34 
35 // forward declarations
36 #ifndef DOXYGEN
37 template <typename number>
38 class Vector;
39 template <typename number>
40 class BlockVector;
41 template <typename number>
42 class FullMatrix;
43 template <typename number>
44 class SparseMatrix;
45 #endif
46 
58 template <typename number>
59 class LAPACKFullMatrix : public TransposeTable<number>
60 {
61 public:
65  using size_type = std::make_unsigned<types::blas_int>::type;
66 
76  explicit LAPACKFullMatrix(const size_type size = 0);
77 
78 
83  LAPACKFullMatrix(const size_type rows, const size_type cols);
84 
85 
96 
102 
109  template <typename number2>
112 
119  template <typename number2>
122 
129  operator=(const number d);
130 
135  operator*=(const number factor);
136 
141  operator/=(const number factor);
142 
153  void
154  set(const size_type i, const size_type j, const number value);
155 
160  void
161  add(const number a, const LAPACKFullMatrix<number> &B);
162 
175  void
176  rank1_update(const number a, const Vector<number> &v);
177 
192  void
193  apply_givens_rotation(const std::array<number, 3> &csr,
194  const size_type i,
195  const size_type k,
196  const bool left = true);
197 
204  template <typename MatrixType>
205  void
206  copy_from(const MatrixType &);
207 
213  void
214  reinit(const size_type size);
215 
238  void
240 
260  void
261  remove_row_and_column(const size_type row, const size_type col);
262 
268  void
269  reinit(const size_type rows, const size_type cols);
270 
274  void
276 
282  size_type
283  m() const;
284 
290  size_type
291  n() const;
292 
306  template <typename MatrixType>
307  void
308  fill(const MatrixType &src,
309  const size_type dst_offset_i = 0,
310  const size_type dst_offset_j = 0,
311  const size_type src_offset_i = 0,
312  const size_type src_offset_j = 0,
313  const number factor = 1.,
314  const bool transpose = false);
315 
316 
344  template <typename number2>
345  void
346  vmult(Vector<number2> & w,
347  const Vector<number2> &v,
348  const bool adding = false) const;
349 
353  void
355  const Vector<number> &v,
356  const bool adding = false) const;
357 
364  template <typename number2>
365  void
366  vmult_add(Vector<number2> &w, const Vector<number2> &v) const;
367 
371  void
372  vmult_add(Vector<number> &w, const Vector<number> &v) const;
373 
385  template <typename number2>
386  void
387  Tvmult(Vector<number2> & w,
388  const Vector<number2> &v,
389  const bool adding = false) const;
390 
394  void
396  const Vector<number> &v,
397  const bool adding = false) const;
398 
405  template <typename number2>
406  void
407  Tvmult_add(Vector<number2> &w, const Vector<number2> &v) const;
408 
412  void
413  Tvmult_add(Vector<number> &w, const Vector<number> &v) const;
414 
415 
430  void
432  const LAPACKFullMatrix<number> &B,
433  const bool adding = false) const;
434 
439  void
441  const LAPACKFullMatrix<number> &B,
442  const bool adding = false) const;
443 
458  void
460  const LAPACKFullMatrix<number> &B,
461  const bool adding = false) const;
462 
467  void
469  const LAPACKFullMatrix<number> &B,
470  const bool adding = false) const;
471 
487  void
489  const LAPACKFullMatrix<number> &B,
490  const Vector<number> & V,
491  const bool adding = false) const;
492 
507  void
509  const LAPACKFullMatrix<number> &B,
510  const bool adding = false) const;
511 
516  void
518  const LAPACKFullMatrix<number> &B,
519  const bool adding = false) const;
520 
536  void
538  const LAPACKFullMatrix<number> &B,
539  const bool adding = false) const;
540 
545  void
547  const LAPACKFullMatrix<number> &B,
548  const bool adding = false) const;
549 
559  void
561 
567  void
568  scale_rows(const Vector<number> &V);
569 
573  void
575 
582  void
584 
604  number
605  reciprocal_condition_number(const number l1_norm) const;
606 
614  number
616 
622  number
623  determinant() const;
624 
628  number
629  l1_norm() const;
630 
634  number
635  linfty_norm() const;
636 
640  number
641  frobenius_norm() const;
642 
647  number
648  trace() const;
649 
655  void
656  invert();
657 
666  void
667  solve(Vector<number> &v, const bool transposed = false) const;
668 
673  void
674  solve(LAPACKFullMatrix<number> &B, const bool transposed = false) const;
675 
694  void
695  compute_eigenvalues(const bool right_eigenvectors = false,
696  const bool left_eigenvectors = false);
697 
717  void
719  const number upper_bound,
720  const number abs_accuracy,
723 
750  void
753  const number lower_bound,
754  const number upper_bound,
755  const number abs_accuracy,
757  std::vector<Vector<number>> &eigenvectors,
758  const types::blas_int itype = 1);
759 
775  void
778  std::vector<Vector<number>> &eigenvectors,
779  const types::blas_int itype = 1);
780 
800  void
801  compute_svd();
802 
822  void
823  compute_inverse_svd(const double threshold = 0.);
824 
829  void
830  compute_inverse_svd_with_kernel(const unsigned int kernel_size);
831 
835  std::complex<number>
836  eigenvalue(const size_type i) const;
837 
842  number
843  singular_value(const size_type i) const;
844 
849  inline const LAPACKFullMatrix<number> &
850  get_svd_u() const;
851 
856  inline const LAPACKFullMatrix<number> &
857  get_svd_vt() const;
858 
887  void
888  print_formatted(std::ostream & out,
889  const unsigned int precision = 3,
890  const bool scientific = true,
891  const unsigned int width = 0,
892  const char * zero_string = " ",
893  const double denominator = 1.,
894  const double threshold = 0.) const;
895 
896 private:
900  number
901  norm(const char type) const;
902 
908 
914 
918  mutable std::vector<number> work;
919 
923  mutable std::vector<types::blas_int> iwork;
924 
931  std::vector<types::blas_int> ipiv;
932 
936  std::vector<number> inv_work;
937 
942  std::vector<typename numbers::NumberTraits<number>::real_type> wr;
943 
948  std::vector<number> wi;
949 
953  std::vector<number> vl;
954 
958  std::vector<number> vr;
959 
964  std::unique_ptr<LAPACKFullMatrix<number>> svd_u;
965 
970  std::unique_ptr<LAPACKFullMatrix<number>> svd_vt;
971 
975  mutable std::mutex mutex;
976 };
977 
978 
979 
985 template <typename number>
987 {
988 public:
989  void
990  initialize(const LAPACKFullMatrix<number> &);
991  void
992  initialize(const LAPACKFullMatrix<number> &, VectorMemory<Vector<number>> &);
993  void
994  vmult(Vector<number> &, const Vector<number> &) const;
995  void
996  Tvmult(Vector<number> &, const Vector<number> &) const;
997  void
999  void
1001 
1002 private:
1005 };
1006 
1007 /*---------------------- Inline functions -----------------------------------*/
1008 
1009 template <typename number>
1010 inline void
1012  const size_type j,
1013  const number value)
1014 {
1015  (*this)(i, j) = value;
1016 }
1017 
1018 
1019 template <typename number>
1022 {
1023  return static_cast<size_type>(this->n_rows());
1024 }
1025 
1026 template <typename number>
1029 {
1030  return static_cast<size_type>(this->n_cols());
1031 }
1032 
1033 template <typename number>
1034 template <typename MatrixType>
1035 inline void
1037 {
1038  this->reinit(M.m(), M.n());
1039 
1040  // loop over the elements of the argument matrix row by row, as suggested
1041  // in the documentation of the sparse matrix iterator class, and
1042  // copy them into the current object
1043  for (size_type row = 0; row < M.m(); ++row)
1044  {
1045  const typename MatrixType::const_iterator end_row = M.end(row);
1046  for (typename MatrixType::const_iterator entry = M.begin(row);
1047  entry != end_row;
1048  ++entry)
1049  this->el(row, entry->column()) = entry->value();
1050  }
1051 
1053 }
1054 
1055 
1056 
1057 template <typename number>
1058 template <typename MatrixType>
1059 inline void
1061  const size_type dst_offset_i,
1062  const size_type dst_offset_j,
1063  const size_type src_offset_i,
1064  const size_type src_offset_j,
1065  const number factor,
1066  const bool transpose)
1067 {
1068  // loop over the elements of the argument matrix row by row, as suggested
1069  // in the documentation of the sparse matrix iterator class
1070  for (size_type row = src_offset_i; row < M.m(); ++row)
1071  {
1072  const typename MatrixType::const_iterator end_row = M.end(row);
1073  for (typename MatrixType::const_iterator entry = M.begin(row);
1074  entry != end_row;
1075  ++entry)
1076  {
1077  const size_type i = transpose ? entry->column() : row;
1078  const size_type j = transpose ? row : entry->column();
1079 
1080  const size_type dst_i = dst_offset_i + i - src_offset_i;
1081  const size_type dst_j = dst_offset_j + j - src_offset_j;
1082  if (dst_i < this->n_rows() && dst_j < this->n_cols())
1083  (*this)(dst_i, dst_j) = factor * entry->value();
1084  }
1085  }
1086 
1088 }
1089 
1090 
1091 template <typename number>
1092 template <typename number2>
1093 void
1095  const Vector<number2> &,
1096  const bool) const
1097 {
1098  Assert(false,
1099  ExcMessage("LAPACKFullMatrix<number>::vmult must be called with a "
1100  "matching Vector<double> vector type."));
1101 }
1102 
1103 
1104 template <typename number>
1105 template <typename number2>
1106 void
1108  const Vector<number2> &) const
1109 {
1110  Assert(false,
1111  ExcMessage("LAPACKFullMatrix<number>::vmult_add must be called with a "
1112  "matching Vector<double> vector type."));
1113 }
1114 
1115 
1116 template <typename number>
1117 template <typename number2>
1118 void
1120  const Vector<number2> &,
1121  const bool) const
1122 {
1123  Assert(false,
1124  ExcMessage("LAPACKFullMatrix<number>::Tvmult must be called with a "
1125  "matching Vector<double> vector type."));
1126 }
1127 
1128 
1129 template <typename number>
1130 template <typename number2>
1131 void
1133  const Vector<number2> &) const
1134 {
1135  Assert(false,
1136  ExcMessage(
1137  "LAPACKFullMatrix<number>::Tvmult_add must be called with a "
1138  "matching Vector<double> vector type."));
1139 }
1140 
1141 
1142 template <typename number>
1143 inline std::complex<number>
1145 {
1147  Assert(wr.size() == this->n_rows(), ExcInternalError());
1148  Assert(wi.size() == this->n_rows(), ExcInternalError());
1149  AssertIndexRange(i, this->n_rows());
1150 
1152  return std::complex<number>(wi[i]);
1153  else
1154  return std::complex<number>(wr[i], wi[i]);
1155 }
1156 
1157 
1158 template <typename number>
1159 inline number
1161 {
1164  AssertIndexRange(i, wr.size());
1165 
1166  return wr[i];
1167 }
1168 
1169 
1170 template <typename number>
1171 inline const LAPACKFullMatrix<number> &
1173 {
1176 
1177  return *svd_u;
1178 }
1179 
1180 
1181 template <typename number>
1182 inline const LAPACKFullMatrix<number> &
1184 {
1187 
1188  return *svd_vt;
1189 }
1190 
1191 
1192 
1194 
1195 #endif
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1075
LAPACKFullMatrix(const size_type size=0)
std::vector< number > work
void rank1_update(const number a, const Vector< number > &v)
void TmTmult(LAPACKFullMatrix< number > &C, const LAPACKFullMatrix< number > &B, const bool adding=false) const
std::unique_ptr< LAPACKFullMatrix< number > > svd_vt
LAPACKFullMatrix< number > & operator/=(const number factor)
Contents is actually a matrix.
std::array< std::pair< Number, Tensor< 1, dim, Number > >, std::integral_constant< int, dim >::value > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
void Tvmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
LAPACKSupport::State state
void remove_row_and_column(const size_type row, const size_type col)
std::vector< types::blas_int > ipiv
size_type n() const
const TableIndices< N > & size() const
size_type n_cols() const
std::complex< number > eigenvalue(const size_type i) const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1648
size_type m() const
size_type n_rows() const
void apply_givens_rotation(const std::array< number, 3 > &csr, const size_type i, const size_type k, const bool left=true)
static const char V
std::vector< number > vr
void compute_eigenvalues_symmetric(const number lower_bound, const number upper_bound, const number abs_accuracy, Vector< number > &eigenvalues, FullMatrix< number > &eigenvectors)
void reinit(const size_type size)
number norm(const char type) const
const LAPACKFullMatrix< number > & get_svd_u() const
void Tmmult(LAPACKFullMatrix< number > &C, const LAPACKFullMatrix< number > &B, const bool adding=false) const
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
std::vector< types::blas_int > iwork
static ::ExceptionBase & ExcState(State arg1)
void set(const size_type i, const size_type j, const number value)
static ::ExceptionBase & ExcInvalidState()
static ::ExceptionBase & ExcMessage(std::string arg1)
void mTmult(LAPACKFullMatrix< number > &C, const LAPACKFullMatrix< number > &B, const bool adding=false) const
#define Assert(cond, exc)
Definition: exceptions.h:1423
void compute_inverse_svd_with_kernel(const unsigned int kernel_size)
void copy_from(const MatrixType &)
std::make_unsigned< types::blas_int >::type size_type
void add(const number a, const LAPACKFullMatrix< number > &B)
std::vector< number > inv_work
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:369
void transpose(LAPACKFullMatrix< number > &B) const
const LAPACKFullMatrix< number > & get_svd_vt() const
Matrix is the inverse of a singular value decomposition.
std::vector< number > vl
std::vector< number > wi
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
std::unique_ptr< LAPACKFullMatrix< number > > svd_u
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void print_formatted(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const unsigned int width=0, const char *zero_string=" ", const double denominator=1., const double threshold=0.) const
void grow_or_shrink(const size_type size)
void scale_rows(const Vector< number > &V)
number singular_value(const size_type i) const
LAPACKFullMatrix< number > & operator*=(const number factor)
void compute_inverse_svd(const double threshold=0.)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SmartPointer< const LAPACKFullMatrix< number >, PreconditionLU< number > > matrix
number linfty_norm() const
LAPACKFullMatrix< number > & operator=(const LAPACKFullMatrix< number > &)
void fill(const MatrixType &src, const size_type dst_offset_i=0, const size_type dst_offset_j=0, const size_type src_offset_i=0, const size_type src_offset_j=0, const number factor=1., const bool transpose=false)
number determinant() const
typename AlignedVector< Number >::size_type size_type
Definition: table.h:409
number frobenius_norm() const
void Tvmult_add(Vector< number2 > &w, const Vector< number2 > &v) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:368
std::vector< typename numbers::NumberTraits< number >::real_type > wr
Matrix contains singular value decomposition,.
void set_property(const LAPACKSupport::Property property)
void solve(Vector< number > &v, const bool transposed=false) const
reference el(const size_type i, const size_type j)
void compute_generalized_eigenvalues_symmetric(LAPACKFullMatrix< number > &B, const number lower_bound, const number upper_bound, const number abs_accuracy, Vector< number > &eigenvalues, std::vector< Vector< number >> &eigenvectors, const types::blas_int itype=1)
void compute_eigenvalues(const bool right_eigenvectors=false, const bool left_eigenvectors=false)
number l1_norm() const
SmartPointer< VectorMemory< Vector< number > >, PreconditionLU< number > > mem
number reciprocal_condition_number() const
Eigenvalue vector is filled.
LAPACKSupport::Property property
void vmult_add(Vector< number2 > &w, const Vector< number2 > &v) const
void mmult(LAPACKFullMatrix< number > &C, const LAPACKFullMatrix< number > &B, const bool adding=false) const
static ::ExceptionBase & ExcInternalError()