Reference documentation for deal.II version GIT 112f7bbc69 2023-05-28 21:25:02+00:00
\(\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 
22 #include <deal.II/base/mutex.h>
24 #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 
82  LAPACKFullMatrix(const size_type rows, const size_type cols);
83 
94 
100 
107  template <typename number2>
110 
117  template <typename number2>
120 
127  operator=(const number d);
128 
133  operator*=(const number factor);
134 
139  operator/=(const number factor);
140 
151  void
152  set(const size_type i, const size_type j, const number value);
153 
158  void
159  add(const number a, const LAPACKFullMatrix<number> &B);
160 
173  void
174  rank1_update(const number a, const Vector<number> &v);
175 
190  void
191  apply_givens_rotation(const std::array<number, 3> &csr,
192  const size_type i,
193  const size_type k,
194  const bool left = true);
195 
202  template <typename MatrixType>
203  void
204  copy_from(const MatrixType &);
205 
211  void
212  reinit(const size_type size);
213 
236  void
238 
258  void
259  remove_row_and_column(const size_type row, const size_type col);
260 
266  void
267  reinit(const size_type rows, const size_type cols);
268 
272  void
274 
280  size_type
281  m() const;
282 
288  size_type
289  n() const;
290 
304  template <typename MatrixType>
305  void
306  fill(const MatrixType &src,
307  const size_type dst_offset_i = 0,
308  const size_type dst_offset_j = 0,
309  const size_type src_offset_i = 0,
310  const size_type src_offset_j = 0,
311  const number factor = 1.,
312  const bool transpose = false);
313 
314 
342  template <typename number2>
343  void
345  const Vector<number2> &v,
346  const bool adding = false) const;
347 
351  void
353  const Vector<number> &v,
354  const bool adding = false) const;
355 
362  template <typename number2>
363  void
365 
369  void
370  vmult_add(Vector<number> &w, const Vector<number> &v) const;
371 
383  template <typename number2>
384  void
386  const Vector<number2> &v,
387  const bool adding = false) const;
388 
392  void
394  const Vector<number> &v,
395  const bool adding = false) const;
396 
403  template <typename number2>
404  void
406 
410  void
411  Tvmult_add(Vector<number> &w, const Vector<number> &v) const;
412 
413 
428  void
430  const LAPACKFullMatrix<number> &B,
431  const bool adding = false) const;
432 
437  void
439  const LAPACKFullMatrix<number> &B,
440  const bool adding = false) const;
441 
456  void
458  const LAPACKFullMatrix<number> &B,
459  const bool adding = false) const;
460 
465  void
467  const LAPACKFullMatrix<number> &B,
468  const bool adding = false) const;
469 
486  void
488  const LAPACKFullMatrix<number> &B,
489  const Vector<number> & V,
490  const bool adding = false) const;
491 
506  void
508  const LAPACKFullMatrix<number> &B,
509  const bool adding = false) const;
510 
515  void
517  const LAPACKFullMatrix<number> &B,
518  const bool adding = false) const;
519 
535  void
537  const LAPACKFullMatrix<number> &B,
538  const bool adding = false) const;
539 
544  void
546  const LAPACKFullMatrix<number> &B,
547  const bool adding = false) const;
548 
558  void
560 
566  void
567  scale_rows(const Vector<number> &V);
568 
572  void
574 
581  void
583 
603  number
604  reciprocal_condition_number(const number l1_norm) const;
605 
613  number
615 
621  number
622  determinant() const;
623 
627  number
628  l1_norm() const;
629 
633  number
634  linfty_norm() const;
635 
639  number
640  frobenius_norm() const;
641 
646  number
647  trace() const;
648 
654  void
655  invert();
656 
665  void
666  solve(Vector<number> &v, const bool transposed = false) const;
667 
672  void
673  solve(LAPACKFullMatrix<number> &B, const bool transposed = false) const;
674 
693  void
694  compute_eigenvalues(const bool right_eigenvectors = false,
695  const bool left_eigenvectors = false);
696 
716  void
718  const number upper_bound,
719  const number abs_accuracy,
722 
749  void
752  const number lower_bound,
753  const number upper_bound,
754  const number abs_accuracy,
756  std::vector<Vector<number>> &eigenvectors,
757  const types::blas_int itype = 1);
758 
774  void
777  std::vector<Vector<number>> &eigenvectors,
778  const types::blas_int itype = 1);
779 
799  void
800  compute_svd();
801 
821  void
822  compute_inverse_svd(const double threshold = 0.);
823 
828  void
829  compute_inverse_svd_with_kernel(const unsigned int kernel_size);
830 
837  std::complex<typename numbers::NumberTraits<number>::real_type>
838  eigenvalue(const size_type i) const;
839 
850  get_right_eigenvectors() const;
851 
858  get_left_eigenvectors() const;
859 
864  number
865  singular_value(const size_type i) const;
866 
871  inline const LAPACKFullMatrix<number> &
872  get_svd_u() const;
873 
878  inline const LAPACKFullMatrix<number> &
879  get_svd_vt() const;
880 
909  void
910  print_formatted(std::ostream & out,
911  const unsigned int precision = 3,
912  const bool scientific = true,
913  const unsigned int width = 0,
914  const char * zero_string = " ",
915  const double denominator = 1.,
916  const double threshold = 0.) const;
917 
918 private:
922  number
923  norm(const char type) const;
924 
930 
936 
940  mutable std::vector<number> work;
941 
945  mutable std::vector<types::blas_int> iwork;
946 
953  std::vector<types::blas_int> ipiv;
954 
958  std::vector<number> inv_work;
959 
964  std::vector<typename numbers::NumberTraits<number>::real_type> wr;
965 
970  std::vector<number> wi;
971 
975  std::vector<number> vl;
976 
980  std::vector<number> vr;
981 
986  std::unique_ptr<LAPACKFullMatrix<number>> svd_u;
987 
992  std::unique_ptr<LAPACKFullMatrix<number>> svd_vt;
993 
998 };
999 
1000 
1001 
1007 template <typename number>
1009 {
1010 public:
1011  void
1013  void
1015  void
1016  vmult(Vector<number> &, const Vector<number> &) const;
1017  void
1018  Tvmult(Vector<number> &, const Vector<number> &) const;
1019  void
1020  vmult(BlockVector<number> &, const BlockVector<number> &) const;
1021  void
1023 
1024 private:
1027 };
1028 
1029 /*---------------------- Inline functions -----------------------------------*/
1030 
1031 template <typename number>
1032 inline void
1034  const size_type j,
1035  const number value)
1036 {
1037  (*this)(i, j) = value;
1038 }
1039 
1040 
1041 
1042 template <typename number>
1045 {
1046  return static_cast<size_type>(this->n_rows());
1047 }
1048 
1049 
1050 
1051 template <typename number>
1054 {
1055  return static_cast<size_type>(this->n_cols());
1056 }
1057 
1058 
1059 
1060 template <typename number>
1061 template <typename MatrixType>
1062 inline void
1064 {
1065  this->reinit(M.m(), M.n());
1066 
1067  // loop over the elements of the argument matrix row by row, as suggested
1068  // in the documentation of the sparse matrix iterator class, and
1069  // copy them into the current object
1070  for (size_type row = 0; 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  this->el(row, entry->column()) = entry->value();
1077  }
1078 
1079  state = LAPACKSupport::matrix;
1080 }
1081 
1082 
1083 
1084 template <typename number>
1085 template <typename MatrixType>
1086 inline void
1088  const size_type dst_offset_i,
1089  const size_type dst_offset_j,
1090  const size_type src_offset_i,
1091  const size_type src_offset_j,
1092  const number factor,
1093  const bool transpose)
1094 {
1095  // loop over the elements of the argument matrix row by row, as suggested
1096  // in the documentation of the sparse matrix iterator class
1097  for (size_type row = src_offset_i; row < M.m(); ++row)
1098  {
1099  const typename MatrixType::const_iterator end_row = M.end(row);
1100  for (typename MatrixType::const_iterator entry = M.begin(row);
1101  entry != end_row;
1102  ++entry)
1103  {
1104  const size_type i = transpose ? entry->column() : row;
1105  const size_type j = transpose ? row : entry->column();
1106 
1107  const size_type dst_i = dst_offset_i + i - src_offset_i;
1108  const size_type dst_j = dst_offset_j + j - src_offset_j;
1109  if (dst_i < this->n_rows() && dst_j < this->n_cols())
1110  (*this)(dst_i, dst_j) = factor * entry->value();
1111  }
1112  }
1113 
1114  state = LAPACKSupport::matrix;
1115 }
1116 
1117 
1118 
1119 template <typename number>
1120 template <typename number2>
1121 void
1123  const Vector<number2> &,
1124  const bool) const
1125 {
1126  Assert(false,
1127  ExcMessage("LAPACKFullMatrix<number>::vmult must be called with a "
1128  "matching Vector<double> vector type."));
1129 }
1130 
1131 
1132 
1133 template <typename number>
1134 template <typename number2>
1135 void
1137  const Vector<number2> &) const
1138 {
1139  Assert(false,
1140  ExcMessage("LAPACKFullMatrix<number>::vmult_add must be called with a "
1141  "matching Vector<double> vector type."));
1142 }
1143 
1144 
1145 
1146 template <typename number>
1147 template <typename number2>
1148 void
1150  const Vector<number2> &,
1151  const bool) const
1152 {
1153  Assert(false,
1154  ExcMessage("LAPACKFullMatrix<number>::Tvmult must be called with a "
1155  "matching Vector<double> vector type."));
1156 }
1157 
1158 
1159 
1160 template <typename number>
1161 template <typename number2>
1162 void
1164  const Vector<number2> &) const
1165 {
1166  Assert(false,
1167  ExcMessage("LAPACKFullMatrix<number>::Tvmult_add must be called "
1168  "with a matching Vector<double> vector type."));
1169 }
1170 
1171 
1172 
1173 namespace internal
1174 {
1175  namespace LAPACKFullMatrixImplementation
1176  {
1177  template <typename RealNumber>
1178  std::complex<RealNumber>
1179  pack_complex(const RealNumber &real_part, const RealNumber &imaginary_part)
1180  {
1181  return std::complex<RealNumber>(real_part, imaginary_part);
1182  }
1183 
1184  // The eigenvalues in LAPACKFullMatrix with complex-valued matrices are
1185  // contained in the 'wi' array, ignoring the 'wr' array.
1186  template <typename Number>
1187  std::complex<Number>
1188  pack_complex(const Number &, const std::complex<Number> &complex_number)
1189  {
1190  return complex_number;
1191  }
1192  } // namespace LAPACKFullMatrixImplementation
1193 } // namespace internal
1194 
1195 
1196 
1197 template <typename number>
1198 inline std::complex<typename numbers::NumberTraits<number>::real_type>
1200 {
1202  Assert(wr.size() == this->n_rows(), ExcInternalError());
1203  Assert(wi.size() == this->n_rows(), ExcInternalError());
1204  AssertIndexRange(i, this->n_rows());
1205 
1207 }
1208 
1209 
1210 
1211 template <typename number>
1212 inline number
1214 {
1216  LAPACKSupport::ExcState(state));
1217  AssertIndexRange(i, wr.size());
1218 
1219  return wr[i];
1220 }
1221 
1222 
1223 
1224 template <typename number>
1225 inline const LAPACKFullMatrix<number> &
1227 {
1229  LAPACKSupport::ExcState(state));
1230 
1231  return *svd_u;
1232 }
1233 
1234 
1235 
1236 template <typename number>
1237 inline const LAPACKFullMatrix<number> &
1239 {
1241  LAPACKSupport::ExcState(state));
1242 
1243  return *svd_vt;
1244 }
1245 
1246 
1247 
1249 
1250 #endif
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
LAPACKFullMatrix< number > & operator*=(const number factor)
number reciprocal_condition_number() const
void Tmmult(LAPACKFullMatrix< number > &C, const LAPACKFullMatrix< number > &B, const bool adding=false) const
void copy_from(const MatrixType &)
void scale_rows(const Vector< number > &V)
FullMatrix< std::complex< typename numbers::NumberTraits< number >::real_type > > get_right_eigenvectors() const
void add(const number a, const LAPACKFullMatrix< number > &B)
void Tvmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
void transpose(LAPACKFullMatrix< number > &B) const
void mTmult(LAPACKFullMatrix< number > &C, const LAPACKFullMatrix< number > &B, const bool adding=false) const
std::vector< typename numbers::NumberTraits< number >::real_type > wr
void compute_eigenvalues_symmetric(const number lower_bound, const number upper_bound, const number abs_accuracy, Vector< number > &eigenvalues, FullMatrix< number > &eigenvectors)
const LAPACKFullMatrix< number > & get_svd_u() const
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
void reinit(const size_type size)
std::make_unsigned< types::blas_int >::type size_type
LAPACKFullMatrix< number > & operator=(const LAPACKFullMatrix< number > &)
number l1_norm() const
FullMatrix< std::complex< typename numbers::NumberTraits< number >::real_type > > get_left_eigenvectors() const
std::unique_ptr< LAPACKFullMatrix< number > > svd_vt
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
std::vector< number > work
void grow_or_shrink(const size_type size)
void apply_givens_rotation(const std::array< number, 3 > &csr, const size_type i, const size_type k, const bool left=true)
void set_property(const LAPACKSupport::Property property)
std::complex< typename numbers::NumberTraits< number >::real_type > eigenvalue(const size_type i) const
number norm(const char type) const
void solve(Vector< number > &v, const bool transposed=false) const
void compute_eigenvalues(const bool right_eigenvectors=false, const bool left_eigenvectors=false)
LAPACKSupport::State state
std::vector< number > inv_work
number frobenius_norm() const
LAPACKFullMatrix(const size_type size=0)
LAPACKSupport::Property property
std::vector< number > wi
size_type m() const
number singular_value(const size_type i) const
void set(const size_type i, const size_type j, const number value)
void compute_inverse_svd(const double threshold=0.)
void vmult_add(Vector< number2 > &w, const Vector< number2 > &v) const
size_type n() const
number linfty_norm() const
void TmTmult(LAPACKFullMatrix< number > &C, const LAPACKFullMatrix< number > &B, const bool adding=false) const
const LAPACKFullMatrix< number > & get_svd_vt() const
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)
std::unique_ptr< LAPACKFullMatrix< number > > svd_u
std::vector< number > vr
void compute_inverse_svd_with_kernel(const unsigned int kernel_size)
void Tvmult_add(Vector< number2 > &w, const Vector< number2 > &v) const
Threads::Mutex mutex
std::vector< types::blas_int > iwork
void rank1_update(const number a, const Vector< number > &v)
std::vector< types::blas_int > ipiv
void remove_row_and_column(const size_type row, const size_type col)
LAPACKFullMatrix< number > & operator/=(const number factor)
number determinant() const
std::vector< number > vl
void mmult(LAPACKFullMatrix< number > &C, const LAPACKFullMatrix< number > &B, const bool adding=false) const
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)
void vmult(Vector< number > &, const Vector< number > &) const
SmartPointer< VectorMemory< Vector< number > >, PreconditionLU< number > > mem
void initialize(const LAPACKFullMatrix< number > &)
SmartPointer< const LAPACKFullMatrix< number >, PreconditionLU< number > > matrix
void Tvmult(Vector< number > &, const Vector< number > &) const
typename AlignedVector< T >::size_type size_type
Definition: table.h:449
const TableIndices< N > & size() const
typename TableBase< 2, T >::size_type size_type
Definition: table.h:1981
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcState(State arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1614
static ::ExceptionBase & ExcInvalidState()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1855
static ::ExceptionBase & ExcMessage(std::string arg1)
@ matrix
Contents is actually a matrix.
@ svd
Matrix contains singular value decomposition,.
@ inverse_svd
Matrix is the inverse of a singular value decomposition.
@ eigenvalues
Eigenvalue vector is filled.
static const char V
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1016
std::complex< RealNumber > pack_complex(const RealNumber &real_part, const RealNumber &imaginary_part)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
int blas_int
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)