Reference documentation for deal.II version GIT c1ddcf411d 2023-11-30 16:30: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 - 2023 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_t<types::blas_int>;
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 
911  void
912  print_formatted(std::ostream &out,
913  const unsigned int precision = 3,
914  const bool scientific = true,
915  const unsigned int width = 0,
916  const char *zero_string = " ",
917  const double denominator = 1.,
918  const double threshold = 0.,
919  const char *separator = " ") const;
920 
921 private:
925  number
926  norm(const char type) const;
927 
933 
939 
943  mutable std::vector<number> work;
944 
948  mutable std::vector<types::blas_int> iwork;
949 
956  std::vector<types::blas_int> ipiv;
957 
961  std::vector<number> inv_work;
962 
967  std::vector<typename numbers::NumberTraits<number>::real_type> wr;
968 
973  std::vector<number> wi;
974 
978  std::vector<number> vl;
979 
983  std::vector<number> vr;
984 
989  std::unique_ptr<LAPACKFullMatrix<number>> svd_u;
990 
995  std::unique_ptr<LAPACKFullMatrix<number>> svd_vt;
996 
1001 };
1002 
1003 
1004 
1010 template <typename number>
1012 {
1013 public:
1014  void
1016  void
1018  void
1019  vmult(Vector<number> &, const Vector<number> &) const;
1020  void
1021  Tvmult(Vector<number> &, const Vector<number> &) const;
1022  void
1023  vmult(BlockVector<number> &, const BlockVector<number> &) const;
1024  void
1026 
1027 private:
1030 };
1031 
1032 /*---------------------- Inline functions -----------------------------------*/
1033 
1034 template <typename number>
1035 inline void
1037  const size_type j,
1038  const number value)
1039 {
1040  (*this)(i, j) = value;
1041 }
1042 
1043 
1044 
1045 template <typename number>
1048 {
1049  return static_cast<size_type>(this->n_rows());
1050 }
1051 
1052 
1053 
1054 template <typename number>
1057 {
1058  return static_cast<size_type>(this->n_cols());
1059 }
1060 
1061 
1062 
1063 template <typename number>
1064 template <typename MatrixType>
1065 inline void
1067 {
1068  this->reinit(M.m(), M.n());
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, and
1072  // copy them into the current object
1073  for (size_type row = 0; row < M.m(); ++row)
1074  {
1075  const typename MatrixType::const_iterator end_row = M.end(row);
1076  for (typename MatrixType::const_iterator entry = M.begin(row);
1077  entry != end_row;
1078  ++entry)
1079  this->el(row, entry->column()) = entry->value();
1080  }
1081 
1082  state = LAPACKSupport::matrix;
1083 }
1084 
1085 
1086 
1087 template <typename number>
1088 template <typename MatrixType>
1089 inline void
1091  const size_type dst_offset_i,
1092  const size_type dst_offset_j,
1093  const size_type src_offset_i,
1094  const size_type src_offset_j,
1095  const number factor,
1096  const bool transpose)
1097 {
1098  // loop over the elements of the argument matrix row by row, as suggested
1099  // in the documentation of the sparse matrix iterator class
1100  for (size_type row = src_offset_i; row < M.m(); ++row)
1101  {
1102  const typename MatrixType::const_iterator end_row = M.end(row);
1103  for (typename MatrixType::const_iterator entry = M.begin(row);
1104  entry != end_row;
1105  ++entry)
1106  {
1107  const size_type i = transpose ? entry->column() : row;
1108  const size_type j = transpose ? row : entry->column();
1109 
1110  const size_type dst_i = dst_offset_i + i - src_offset_i;
1111  const size_type dst_j = dst_offset_j + j - src_offset_j;
1112  if (dst_i < this->n_rows() && dst_j < this->n_cols())
1113  (*this)(dst_i, dst_j) = factor * entry->value();
1114  }
1115  }
1116 
1117  state = LAPACKSupport::matrix;
1118 }
1119 
1120 
1121 
1122 template <typename number>
1123 template <typename number2>
1124 void
1126  const Vector<number2> &,
1127  const bool) const
1128 {
1129  Assert(false,
1130  ExcMessage("LAPACKFullMatrix<number>::vmult must be called with a "
1131  "matching Vector<double> vector type."));
1132 }
1133 
1134 
1135 
1136 template <typename number>
1137 template <typename number2>
1138 void
1140  const Vector<number2> &) const
1141 {
1142  Assert(false,
1143  ExcMessage("LAPACKFullMatrix<number>::vmult_add must be called with a "
1144  "matching Vector<double> vector type."));
1145 }
1146 
1147 
1148 
1149 template <typename number>
1150 template <typename number2>
1151 void
1153  const Vector<number2> &,
1154  const bool) const
1155 {
1156  Assert(false,
1157  ExcMessage("LAPACKFullMatrix<number>::Tvmult must be called with a "
1158  "matching Vector<double> vector type."));
1159 }
1160 
1161 
1162 
1163 template <typename number>
1164 template <typename number2>
1165 void
1167  const Vector<number2> &) const
1168 {
1169  Assert(false,
1170  ExcMessage("LAPACKFullMatrix<number>::Tvmult_add must be called "
1171  "with a matching Vector<double> vector type."));
1172 }
1173 
1174 
1175 
1176 namespace internal
1177 {
1178  namespace LAPACKFullMatrixImplementation
1179  {
1180  template <typename RealNumber>
1181  std::complex<RealNumber>
1182  pack_complex(const RealNumber &real_part, const RealNumber &imaginary_part)
1183  {
1184  return std::complex<RealNumber>(real_part, imaginary_part);
1185  }
1186 
1187  // The eigenvalues in LAPACKFullMatrix with complex-valued matrices are
1188  // contained in the 'wi' array, ignoring the 'wr' array.
1189  template <typename Number>
1190  std::complex<Number>
1191  pack_complex(const Number &, const std::complex<Number> &complex_number)
1192  {
1193  return complex_number;
1194  }
1195  } // namespace LAPACKFullMatrixImplementation
1196 } // namespace internal
1197 
1198 
1199 
1200 template <typename number>
1201 inline std::complex<typename numbers::NumberTraits<number>::real_type>
1203 {
1205  Assert(wr.size() == this->n_rows(), ExcInternalError());
1206  Assert(wi.size() == this->n_rows(), ExcInternalError());
1207  AssertIndexRange(i, this->n_rows());
1208 
1210 }
1211 
1212 
1213 
1214 template <typename number>
1215 inline number
1217 {
1219  LAPACKSupport::ExcState(state));
1220  AssertIndexRange(i, wr.size());
1221 
1222  return wr[i];
1223 }
1224 
1225 
1226 
1227 template <typename number>
1228 inline const LAPACKFullMatrix<number> &
1230 {
1232  LAPACKSupport::ExcState(state));
1233 
1234  return *svd_u;
1235 }
1236 
1237 
1238 
1239 template <typename number>
1240 inline const LAPACKFullMatrix<number> &
1242 {
1244  LAPACKSupport::ExcState(state));
1245 
1246  return *svd_vt;
1247 }
1248 
1249 
1250 
1252 
1253 #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)
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
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::make_unsigned_t< types::blas_int > size_type
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 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 char *separator=" ") 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:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcState(State arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcInvalidState()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1888
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:1008
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)