Reference documentation for deal.II version 9.2.0
\(\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 
59 template <typename number>
60 class LAPACKFullMatrix : public TransposeTable<number>
61 {
62 public:
66  using size_type = std::make_unsigned<types::blas_int>::type;
67 
77  explicit LAPACKFullMatrix(const size_type size = 0);
78 
79 
84  LAPACKFullMatrix(const size_type rows, const size_type cols);
85 
86 
97 
103 
110  template <typename number2>
113 
120  template <typename number2>
123 
130  operator=(const number d);
131 
136  operator*=(const number factor);
137 
142  operator/=(const number factor);
143 
154  void
155  set(const size_type i, const size_type j, const number value);
156 
161  void
162  add(const number a, const LAPACKFullMatrix<number> &B);
163 
176  void
177  rank1_update(const number a, const Vector<number> &v);
178 
193  void
194  apply_givens_rotation(const std::array<number, 3> &csr,
195  const size_type i,
196  const size_type k,
197  const bool left = true);
198 
205  template <typename MatrixType>
206  void
207  copy_from(const MatrixType &);
208 
214  void
215  reinit(const size_type size);
216 
239  void
241 
261  void
262  remove_row_and_column(const size_type row, const size_type col);
263 
269  void
270  reinit(const size_type rows, const size_type cols);
271 
275  void
277 
283  size_type
284  m() const;
285 
291  size_type
292  n() const;
293 
307  template <typename MatrixType>
308  void
309  fill(const MatrixType &src,
310  const size_type dst_offset_i = 0,
311  const size_type dst_offset_j = 0,
312  const size_type src_offset_i = 0,
313  const size_type src_offset_j = 0,
314  const number factor = 1.,
315  const bool transpose = false);
316 
317 
345  template <typename number2>
346  void
347  vmult(Vector<number2> & w,
348  const Vector<number2> &v,
349  const bool adding = false) const;
350 
354  void
356  const Vector<number> &v,
357  const bool adding = false) const;
358 
365  template <typename number2>
366  void
367  vmult_add(Vector<number2> &w, const Vector<number2> &v) const;
368 
372  void
373  vmult_add(Vector<number> &w, const Vector<number> &v) const;
374 
386  template <typename number2>
387  void
388  Tvmult(Vector<number2> & w,
389  const Vector<number2> &v,
390  const bool adding = false) const;
391 
395  void
397  const Vector<number> &v,
398  const bool adding = false) const;
399 
406  template <typename number2>
407  void
408  Tvmult_add(Vector<number2> &w, const Vector<number2> &v) const;
409 
413  void
414  Tvmult_add(Vector<number> &w, const Vector<number> &v) const;
415 
416 
431  void
433  const LAPACKFullMatrix<number> &B,
434  const bool adding = false) const;
435 
440  void
442  const LAPACKFullMatrix<number> &B,
443  const bool adding = false) const;
444 
459  void
461  const LAPACKFullMatrix<number> &B,
462  const bool adding = false) const;
463 
468  void
470  const LAPACKFullMatrix<number> &B,
471  const bool adding = false) const;
472 
488  void
490  const LAPACKFullMatrix<number> &B,
491  const Vector<number> & V,
492  const bool adding = false) const;
493 
508  void
510  const LAPACKFullMatrix<number> &B,
511  const bool adding = false) const;
512 
517  void
519  const LAPACKFullMatrix<number> &B,
520  const bool adding = false) const;
521 
537  void
539  const LAPACKFullMatrix<number> &B,
540  const bool adding = false) const;
541 
546  void
548  const LAPACKFullMatrix<number> &B,
549  const bool adding = false) const;
550 
560  void
562 
568  void
569  scale_rows(const Vector<number> &V);
570 
574  void
576 
583  void
585 
605  number
606  reciprocal_condition_number(const number l1_norm) const;
607 
615  number
617 
623  number
624  determinant() const;
625 
629  number
630  l1_norm() const;
631 
635  number
636  linfty_norm() const;
637 
641  number
642  frobenius_norm() const;
643 
648  number
649  trace() const;
650 
656  void
657  invert();
658 
667  void
668  solve(Vector<number> &v, const bool transposed = false) const;
669 
674  void
675  solve(LAPACKFullMatrix<number> &B, const bool transposed = false) const;
676 
695  void
696  compute_eigenvalues(const bool right_eigenvectors = false,
697  const bool left_eigenvectors = false);
698 
718  void
720  const number upper_bound,
721  const number abs_accuracy,
724 
751  void
754  const number lower_bound,
755  const number upper_bound,
756  const number abs_accuracy,
758  std::vector<Vector<number>> &eigenvectors,
759  const types::blas_int itype = 1);
760 
776  void
779  std::vector<Vector<number>> &eigenvectors,
780  const types::blas_int itype = 1);
781 
801  void
802  compute_svd();
803 
823  void
824  compute_inverse_svd(const double threshold = 0.);
825 
830  void
831  compute_inverse_svd_with_kernel(const unsigned int kernel_size);
832 
836  std::complex<number>
837  eigenvalue(const size_type i) const;
838 
843  number
844  singular_value(const size_type i) const;
845 
850  inline const LAPACKFullMatrix<number> &
851  get_svd_u() const;
852 
857  inline const LAPACKFullMatrix<number> &
858  get_svd_vt() const;
859 
888  void
889  print_formatted(std::ostream & out,
890  const unsigned int precision = 3,
891  const bool scientific = true,
892  const unsigned int width = 0,
893  const char * zero_string = " ",
894  const double denominator = 1.,
895  const double threshold = 0.) const;
896 
897 private:
901  number
902  norm(const char type) const;
903 
909 
915 
919  mutable std::vector<number> work;
920 
924  mutable std::vector<types::blas_int> iwork;
925 
932  std::vector<types::blas_int> ipiv;
933 
937  std::vector<number> inv_work;
938 
943  std::vector<typename numbers::NumberTraits<number>::real_type> wr;
944 
949  std::vector<number> wi;
950 
954  std::vector<number> vl;
955 
959  std::vector<number> vr;
960 
965  std::unique_ptr<LAPACKFullMatrix<number>> svd_u;
966 
971  std::unique_ptr<LAPACKFullMatrix<number>> svd_vt;
972 
976  mutable std::mutex mutex;
977 };
978 
979 
980 
987 template <typename number>
989 {
990 public:
991  void
992  initialize(const LAPACKFullMatrix<number> &);
993  void
994  initialize(const LAPACKFullMatrix<number> &, VectorMemory<Vector<number>> &);
995  void
996  vmult(Vector<number> &, const Vector<number> &) const;
997  void
998  Tvmult(Vector<number> &, const Vector<number> &) const;
999  void
1000  vmult(BlockVector<number> &, const BlockVector<number> &) const;
1001  void
1003 
1004 private:
1007 };
1008 
1009 /*---------------------- Inline functions -----------------------------------*/
1010 
1011 template <typename number>
1012 inline void
1014  const size_type j,
1015  const number value)
1016 {
1017  (*this)(i, j) = value;
1018 }
1019 
1020 
1021 template <typename number>
1024 {
1025  return static_cast<size_type>(this->n_rows());
1026 }
1027 
1028 template <typename number>
1031 {
1032  return static_cast<size_type>(this->n_cols());
1033 }
1034 
1035 template <typename number>
1036 template <typename MatrixType>
1037 inline void
1039 {
1040  this->reinit(M.m(), M.n());
1041 
1042  // loop over the elements of the argument matrix row by row, as suggested
1043  // in the documentation of the sparse matrix iterator class, and
1044  // copy them into the current object
1045  for (size_type row = 0; row < M.m(); ++row)
1046  {
1047  const typename MatrixType::const_iterator end_row = M.end(row);
1048  for (typename MatrixType::const_iterator entry = M.begin(row);
1049  entry != end_row;
1050  ++entry)
1051  this->el(row, entry->column()) = entry->value();
1052  }
1053 
1055 }
1056 
1057 
1058 
1059 template <typename number>
1060 template <typename MatrixType>
1061 inline void
1063  const size_type dst_offset_i,
1064  const size_type dst_offset_j,
1065  const size_type src_offset_i,
1066  const size_type src_offset_j,
1067  const number factor,
1068  const bool transpose)
1069 {
1070  // loop over the elements of the argument matrix row by row, as suggested
1071  // in the documentation of the sparse matrix iterator class
1072  for (size_type row = src_offset_i; row < M.m(); ++row)
1073  {
1074  const typename MatrixType::const_iterator end_row = M.end(row);
1075  for (typename MatrixType::const_iterator entry = M.begin(row);
1076  entry != end_row;
1077  ++entry)
1078  {
1079  const size_type i = transpose ? entry->column() : row;
1080  const size_type j = transpose ? row : entry->column();
1081 
1082  const size_type dst_i = dst_offset_i + i - src_offset_i;
1083  const size_type dst_j = dst_offset_j + j - src_offset_j;
1084  if (dst_i < this->n_rows() && dst_j < this->n_cols())
1085  (*this)(dst_i, dst_j) = factor * entry->value();
1086  }
1087  }
1088 
1090 }
1091 
1092 
1093 template <typename number>
1094 template <typename number2>
1095 void
1097  const Vector<number2> &,
1098  const bool) const
1099 {
1100  Assert(false,
1101  ExcMessage("LAPACKFullMatrix<number>::vmult must be called with a "
1102  "matching Vector<double> vector type."));
1103 }
1104 
1105 
1106 template <typename number>
1107 template <typename number2>
1108 void
1110  const Vector<number2> &) const
1111 {
1112  Assert(false,
1113  ExcMessage("LAPACKFullMatrix<number>::vmult_add must be called with a "
1114  "matching Vector<double> vector type."));
1115 }
1116 
1117 
1118 template <typename number>
1119 template <typename number2>
1120 void
1122  const Vector<number2> &,
1123  const bool) const
1124 {
1125  Assert(false,
1126  ExcMessage("LAPACKFullMatrix<number>::Tvmult must be called with a "
1127  "matching Vector<double> vector type."));
1128 }
1129 
1130 
1131 template <typename number>
1132 template <typename number2>
1133 void
1135  const Vector<number2> &) const
1136 {
1137  Assert(false,
1138  ExcMessage(
1139  "LAPACKFullMatrix<number>::Tvmult_add must be called with a "
1140  "matching Vector<double> vector type."));
1141 }
1142 
1143 
1144 template <typename number>
1145 inline std::complex<number>
1147 {
1149  Assert(wr.size() == this->n_rows(), ExcInternalError());
1150  Assert(wi.size() == this->n_rows(), ExcInternalError());
1151  AssertIndexRange(i, this->n_rows());
1152 
1154  return std::complex<number>(wi[i]);
1155  else
1156  return std::complex<number>(wr[i], wi[i]);
1157 }
1158 
1159 
1160 template <typename number>
1161 inline number
1163 {
1166  AssertIndexRange(i, wr.size());
1167 
1168  return wr[i];
1169 }
1170 
1171 
1172 template <typename number>
1173 inline const LAPACKFullMatrix<number> &
1175 {
1178 
1179  return *svd_u;
1180 }
1181 
1182 
1183 template <typename number>
1184 inline const LAPACKFullMatrix<number> &
1186 {
1189 
1190  return *svd_vt;
1191 }
1192 
1193 
1194 
1196 
1197 #endif
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1102
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:1649
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:1419
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:359
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:420
number frobenius_norm() const
void Tvmult_add(Vector< number2 > &w, const Vector< number2 > &v) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
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
static const bool value
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()