Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
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 template <typename number>
37 class Vector;
38 template <typename number>
39 class BlockVector;
40 template <typename number>
41 class FullMatrix;
42 template <typename number>
43 class SparseMatrix;
44 
45 
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
354  vmult(Vector<number> & w,
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 
686  DEAL_II_DEPRECATED
687  void
688  apply_lu_factorization(Vector<number> &v, const bool transposed) const;
689 
701  DEAL_II_DEPRECATED
702  void
704  const bool transposed) const;
705 
724  void
725  compute_eigenvalues(const bool right_eigenvectors = false,
726  const bool left_eigenvectors = false);
727 
747  void
748  compute_eigenvalues_symmetric(const number lower_bound,
749  const number upper_bound,
750  const number abs_accuracy,
753 
780  void
783  const number lower_bound,
784  const number upper_bound,
785  const number abs_accuracy,
787  std::vector<Vector<number>> &eigenvectors,
788  const types::blas_int itype = 1);
789 
805  void
808  std::vector<Vector<number>> &eigenvectors,
809  const types::blas_int itype = 1);
810 
819  void
820  compute_svd();
821 
841  void
842  compute_inverse_svd(const double threshold = 0.);
843 
848  void
849  compute_inverse_svd_with_kernel(const unsigned int kernel_size);
850 
854  std::complex<number>
855  eigenvalue(const size_type i) const;
856 
861  number
862  singular_value(const size_type i) const;
863 
892  void
893  print_formatted(std::ostream & out,
894  const unsigned int precision = 3,
895  const bool scientific = true,
896  const unsigned int width = 0,
897  const char * zero_string = " ",
898  const double denominator = 1.,
899  const double threshold = 0.) const;
900 
901 private:
905  number
906  norm(const char type) const;
907 
913 
919 
923  mutable std::vector<number> work;
924 
928  mutable std::vector<types::blas_int> iwork;
929 
936  std::vector<types::blas_int> ipiv;
937 
941  std::vector<number> inv_work;
942 
947  std::vector<typename numbers::NumberTraits<number>::real_type> wr;
948 
953  std::vector<number> wi;
954 
958  std::vector<number> vl;
959 
963  std::vector<number> vr;
964 
969  std::unique_ptr<LAPACKFullMatrix<number>> svd_u;
970 
975  std::unique_ptr<LAPACKFullMatrix<number>> svd_vt;
976 
981 };
982 
983 
984 
991 template <typename number>
993 {
994 public:
995  void
996  initialize(const LAPACKFullMatrix<number> &);
997  void
998  initialize(const LAPACKFullMatrix<number> &, VectorMemory<Vector<number>> &);
999  void
1000  vmult(Vector<number> &, const Vector<number> &) const;
1001  void
1002  Tvmult(Vector<number> &, const Vector<number> &) const;
1003  void
1004  vmult(BlockVector<number> &, const BlockVector<number> &) const;
1005  void
1006  Tvmult(BlockVector<number> &, const BlockVector<number> &) const;
1007 
1008 private:
1011 };
1012 
1013 /*---------------------- Inline functions -----------------------------------*/
1014 
1015 template <typename number>
1016 inline void
1018  const size_type j,
1019  const number value)
1020 {
1021  (*this)(i, j) = value;
1022 }
1023 
1024 
1025 template <typename number>
1028 {
1029  return static_cast<size_type>(this->n_rows());
1030 }
1031 
1032 template <typename number>
1035 {
1036  return static_cast<size_type>(this->n_cols());
1037 }
1038 
1039 template <typename number>
1040 template <typename MatrixType>
1041 inline void
1043 {
1044  this->reinit(M.m(), M.n());
1045 
1046  // loop over the elements of the argument matrix row by row, as suggested
1047  // in the documentation of the sparse matrix iterator class, and
1048  // copy them into the current object
1049  for (size_type row = 0; row < M.m(); ++row)
1050  {
1051  const typename MatrixType::const_iterator end_row = M.end(row);
1052  for (typename MatrixType::const_iterator entry = M.begin(row);
1053  entry != end_row;
1054  ++entry)
1055  this->el(row, entry->column()) = entry->value();
1056  }
1057 
1058  state = LAPACKSupport::matrix;
1059 }
1060 
1061 
1062 
1063 template <typename number>
1064 template <typename MatrixType>
1065 inline void
1067  const size_type dst_offset_i,
1068  const size_type dst_offset_j,
1069  const size_type src_offset_i,
1070  const size_type src_offset_j,
1071  const number factor,
1072  const bool transpose)
1073 {
1074  // loop over the elements of the argument matrix row by row, as suggested
1075  // in the documentation of the sparse matrix iterator class
1076  for (size_type row = src_offset_i; row < M.m(); ++row)
1077  {
1078  const typename MatrixType::const_iterator end_row = M.end(row);
1079  for (typename MatrixType::const_iterator entry = M.begin(row);
1080  entry != end_row;
1081  ++entry)
1082  {
1083  const size_type i = transpose ? entry->column() : row;
1084  const size_type j = transpose ? row : entry->column();
1085 
1086  const size_type dst_i = dst_offset_i + i - src_offset_i;
1087  const size_type dst_j = dst_offset_j + j - src_offset_j;
1088  if (dst_i < this->n_rows() && dst_j < this->n_cols())
1089  (*this)(dst_i, dst_j) = factor * entry->value();
1090  }
1091  }
1092 
1093  state = LAPACKSupport::matrix;
1094 }
1095 
1096 
1097 template <typename number>
1098 template <typename number2>
1099 void
1101  const Vector<number2> &,
1102  const bool) const
1103 {
1104  Assert(false,
1105  ExcMessage("LAPACKFullMatrix<number>::vmult must be called with a "
1106  "matching Vector<double> vector type."));
1107 }
1108 
1109 
1110 template <typename number>
1111 template <typename number2>
1112 void
1114  const Vector<number2> &) const
1115 {
1116  Assert(false,
1117  ExcMessage("LAPACKFullMatrix<number>::vmult_add must be called with a "
1118  "matching Vector<double> vector type."));
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>::Tvmult must be called with a "
1131  "matching Vector<double> vector type."));
1132 }
1133 
1134 
1135 template <typename number>
1136 template <typename number2>
1137 void
1139  const Vector<number2> &) const
1140 {
1141  Assert(false,
1142  ExcMessage(
1143  "LAPACKFullMatrix<number>::Tvmult_add must be called with a "
1144  "matching Vector<double> vector type."));
1145 }
1146 
1147 
1148 template <typename number>
1149 inline std::complex<number>
1151 {
1153  Assert(wr.size() == this->n_rows(), ExcInternalError());
1154  Assert(wi.size() == this->n_rows(), ExcInternalError());
1155  AssertIndexRange(i, this->n_rows());
1156 
1158  return std::complex<number>(wi[i]);
1159  else
1160  return std::complex<number>(wr[i], wi[i]);
1161 }
1162 
1163 
1164 template <typename number>
1165 inline number
1167 {
1169  LAPACKSupport::ExcState(state));
1170  AssertIndexRange(i, wr.size());
1171 
1172  return wr[i];
1173 }
1174 
1175 
1176 
1177 DEAL_II_NAMESPACE_CLOSE
1178 
1179 #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
typename TableBase< 2, T >::size_type size_type
Definition: table.h:1936
void remove_row_and_column(const size_type row, const size_type col)
Threads::Mutex mutex
std::vector< types::blas_int > ipiv
size_type n() const
const TableIndices< N > & size() const
std::complex< number > eigenvalue(const size_type i) const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1637
size_type m() 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
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:1407
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
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
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
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.
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
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()