Reference documentation for deal.II version Git 32d6d15a92 2020-01-28 17:18:58 -0700
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
lapack_full_matrix.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2005 - 2019 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/smartpointer.h>
23 #include <deal.II/base/table.h>
24 #include <deal.II/base/thread_management.h>
25 
26 #include <deal.II/lac/lapack_support.h>
27 #include <deal.II/lac/vector_memory.h>
28 
29 #include <complex>
30 #include <memory>
31 #include <vector>
32 
33 DEAL_II_NAMESPACE_OPEN
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
355  vmult(Vector<number> & w,
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 
687  DEAL_II_DEPRECATED
688  void
689  apply_lu_factorization(Vector<number> &v, const bool transposed) const;
690 
702  DEAL_II_DEPRECATED
703  void
705  const bool transposed) const;
706 
725  void
726  compute_eigenvalues(const bool right_eigenvectors = false,
727  const bool left_eigenvectors = false);
728 
748  void
749  compute_eigenvalues_symmetric(const number lower_bound,
750  const number upper_bound,
751  const number abs_accuracy,
754 
781  void
784  const number lower_bound,
785  const number upper_bound,
786  const number abs_accuracy,
788  std::vector<Vector<number>> &eigenvectors,
789  const types::blas_int itype = 1);
790 
806  void
809  std::vector<Vector<number>> &eigenvectors,
810  const types::blas_int itype = 1);
811 
831  void
832  compute_svd();
833 
853  void
854  compute_inverse_svd(const double threshold = 0.);
855 
860  void
861  compute_inverse_svd_with_kernel(const unsigned int kernel_size);
862 
866  std::complex<number>
867  eigenvalue(const size_type i) const;
868 
873  number
874  singular_value(const size_type i) const;
875 
880  inline const LAPACKFullMatrix<number> &
881  get_svd_u() const;
882 
887  inline const LAPACKFullMatrix<number> &
888  get_svd_vt() const;
889 
918  void
919  print_formatted(std::ostream & out,
920  const unsigned int precision = 3,
921  const bool scientific = true,
922  const unsigned int width = 0,
923  const char * zero_string = " ",
924  const double denominator = 1.,
925  const double threshold = 0.) const;
926 
927 private:
931  number
932  norm(const char type) const;
933 
939 
945 
949  mutable std::vector<number> work;
950 
954  mutable std::vector<types::blas_int> iwork;
955 
962  std::vector<types::blas_int> ipiv;
963 
967  std::vector<number> inv_work;
968 
973  std::vector<typename numbers::NumberTraits<number>::real_type> wr;
974 
979  std::vector<number> wi;
980 
984  std::vector<number> vl;
985 
989  std::vector<number> vr;
990 
995  std::unique_ptr<LAPACKFullMatrix<number>> svd_u;
996 
1001  std::unique_ptr<LAPACKFullMatrix<number>> svd_vt;
1002 
1006  mutable std::mutex mutex;
1007 };
1008 
1009 
1010 
1017 template <typename number>
1019 {
1020 public:
1021  void
1022  initialize(const LAPACKFullMatrix<number> &);
1023  void
1024  initialize(const LAPACKFullMatrix<number> &, VectorMemory<Vector<number>> &);
1025  void
1026  vmult(Vector<number> &, const Vector<number> &) const;
1027  void
1028  Tvmult(Vector<number> &, const Vector<number> &) const;
1029  void
1030  vmult(BlockVector<number> &, const BlockVector<number> &) const;
1031  void
1033 
1034 private:
1037 };
1038 
1039 /*---------------------- Inline functions -----------------------------------*/
1040 
1041 template <typename number>
1042 inline void
1044  const size_type j,
1045  const number value)
1046 {
1047  (*this)(i, j) = value;
1048 }
1049 
1050 
1051 template <typename number>
1054 {
1055  return static_cast<size_type>(this->n_rows());
1056 }
1057 
1058 template <typename number>
1061 {
1062  return static_cast<size_type>(this->n_cols());
1063 }
1064 
1065 template <typename number>
1066 template <typename MatrixType>
1067 inline void
1069 {
1070  this->reinit(M.m(), M.n());
1071 
1072  // loop over the elements of the argument matrix row by row, as suggested
1073  // in the documentation of the sparse matrix iterator class, and
1074  // copy them into the current object
1075  for (size_type row = 0; row < M.m(); ++row)
1076  {
1077  const typename MatrixType::const_iterator end_row = M.end(row);
1078  for (typename MatrixType::const_iterator entry = M.begin(row);
1079  entry != end_row;
1080  ++entry)
1081  this->el(row, entry->column()) = entry->value();
1082  }
1083 
1085 }
1086 
1087 
1088 
1089 template <typename number>
1090 template <typename MatrixType>
1091 inline void
1093  const size_type dst_offset_i,
1094  const size_type dst_offset_j,
1095  const size_type src_offset_i,
1096  const size_type src_offset_j,
1097  const number factor,
1098  const bool transpose)
1099 {
1100  // loop over the elements of the argument matrix row by row, as suggested
1101  // in the documentation of the sparse matrix iterator class
1102  for (size_type row = src_offset_i; row < M.m(); ++row)
1103  {
1104  const typename MatrixType::const_iterator end_row = M.end(row);
1105  for (typename MatrixType::const_iterator entry = M.begin(row);
1106  entry != end_row;
1107  ++entry)
1108  {
1109  const size_type i = transpose ? entry->column() : row;
1110  const size_type j = transpose ? row : entry->column();
1111 
1112  const size_type dst_i = dst_offset_i + i - src_offset_i;
1113  const size_type dst_j = dst_offset_j + j - src_offset_j;
1114  if (dst_i < this->n_rows() && dst_j < this->n_cols())
1115  (*this)(dst_i, dst_j) = factor * entry->value();
1116  }
1117  }
1118 
1120 }
1121 
1122 
1123 template <typename number>
1124 template <typename number2>
1125 void
1127  const Vector<number2> &,
1128  const bool) const
1129 {
1130  Assert(false,
1131  ExcMessage("LAPACKFullMatrix<number>::vmult must be called with a "
1132  "matching Vector<double> vector type."));
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 template <typename number>
1149 template <typename number2>
1150 void
1152  const Vector<number2> &,
1153  const bool) const
1154 {
1155  Assert(false,
1156  ExcMessage("LAPACKFullMatrix<number>::Tvmult must be called with a "
1157  "matching Vector<double> vector type."));
1158 }
1159 
1160 
1161 template <typename number>
1162 template <typename number2>
1163 void
1165  const Vector<number2> &) const
1166 {
1167  Assert(false,
1168  ExcMessage(
1169  "LAPACKFullMatrix<number>::Tvmult_add must be called with a "
1170  "matching Vector<double> vector type."));
1171 }
1172 
1173 
1174 template <typename number>
1175 inline std::complex<number>
1177 {
1179  Assert(wr.size() == this->n_rows(), ExcInternalError());
1180  Assert(wi.size() == this->n_rows(), ExcInternalError());
1181  AssertIndexRange(i, this->n_rows());
1182 
1184  return std::complex<number>(wi[i]);
1185  else
1186  return std::complex<number>(wr[i], wi[i]);
1187 }
1188 
1189 
1190 template <typename number>
1191 inline number
1193 {
1196  AssertIndexRange(i, wr.size());
1197 
1198  return wr[i];
1199 }
1200 
1201 
1202 template <typename number>
1203 inline const LAPACKFullMatrix<number> &
1205 {
1208 
1209  return *svd_u;
1210 }
1211 
1212 
1213 template <typename number>
1214 inline const LAPACKFullMatrix<number> &
1216 {
1219 
1220  return *svd_vt;
1221 }
1222 
1223 
1224 
1225 DEAL_II_NAMESPACE_CLOSE
1226 
1227 #endif
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
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
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:1641
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)
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
LinearAlgebra::distributed::Vector< Number > Vector
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:1411
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
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
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.)
number linfty_norm() const
LAPACKFullMatrix< number > & operator=(const LAPACKFullMatrix< number > &)
int blas_int
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
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)
void apply_lu_factorization(Vector< number > &v, const bool transposed) const
number l1_norm() const
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()