Reference documentation for deal.II version GIT relicensing-224-gc660c0d696 2024-03-28 18:40: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\}}\)
Loading...
Searching...
No Matches
lapack_full_matrix.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2005 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_lapack_full_matrix_h
16#define dealii_lapack_full_matrix_h
17
18
19#include <deal.II/base/config.h>
20
21#include <deal.II/base/mutex.h>
23#include <deal.II/base/table.h>
24
27
28#include <complex>
29#include <memory>
30#include <vector>
31
33
34// forward declarations
35#ifndef DOXYGEN
36template <typename number>
37class Vector;
38template <typename number>
39class BlockVector;
40template <typename number>
41class FullMatrix;
42template <typename number>
43class SparseMatrix;
44#endif
45
57template <typename number>
58class LAPACKFullMatrix : public TransposeTable<number>
59{
60public:
64 using size_type = std::make_unsigned_t<types::blas_int>;
65
75 explicit LAPACKFullMatrix(const size_type size = 0);
76
81 LAPACKFullMatrix(const size_type rows, const size_type cols);
82
93
99
106 template <typename number2>
109
116 template <typename number2>
119
126 operator=(const number d);
127
132 operator*=(const number factor);
133
138 operator/=(const number factor);
139
150 void
151 set(const size_type i, const size_type j, const number value);
152
157 void
158 add(const number a, const LAPACKFullMatrix<number> &B);
159
172 void
173 rank1_update(const number a, const Vector<number> &v);
174
189 void
190 apply_givens_rotation(const std::array<number, 3> &csr,
191 const size_type i,
192 const size_type k,
193 const bool left = true);
194
201 template <typename MatrixType>
202 void
203 copy_from(const MatrixType &);
204
210 void
211 reinit(const size_type size);
212
235 void
237
257 void
258 remove_row_and_column(const size_type row, const size_type col);
259
265 void
266 reinit(const size_type rows, const size_type cols);
267
271 void
273
280 m() const;
281
288 n() const;
289
303 template <typename MatrixType>
304 void
305 fill(const MatrixType &src,
306 const size_type dst_offset_i = 0,
307 const size_type dst_offset_j = 0,
308 const size_type src_offset_i = 0,
309 const size_type src_offset_j = 0,
310 const number factor = 1.,
311 const bool transpose = false);
312
313
341 template <typename number2>
342 void
344 const Vector<number2> &v,
345 const bool adding = false) const;
346
350 void
352 const Vector<number> &v,
353 const bool adding = false) const;
354
361 template <typename number2>
362 void
364
368 void
369 vmult_add(Vector<number> &w, const Vector<number> &v) const;
370
382 template <typename number2>
383 void
385 const Vector<number2> &v,
386 const bool adding = false) const;
387
391 void
393 const Vector<number> &v,
394 const bool adding = false) const;
395
402 template <typename number2>
403 void
405
409 void
410 Tvmult_add(Vector<number> &w, const Vector<number> &v) const;
411
412
427 void
430 const bool adding = false) const;
431
436 void
439 const bool adding = false) const;
440
455 void
458 const bool adding = false) const;
459
464 void
467 const bool adding = false) const;
468
485 void
488 const Vector<number> &V,
489 const bool adding = false) const;
490
505 void
508 const bool adding = false) const;
509
514 void
517 const bool adding = false) const;
518
534 void
537 const bool adding = false) const;
538
543 void
546 const bool adding = false) const;
547
557 void
559
565 void
566 scale_rows(const Vector<number> &V);
567
571 void
573
580 void
582
602 number
603 reciprocal_condition_number(const number l1_norm) const;
604
612 number
614
620 number
621 determinant() const;
622
626 number
627 l1_norm() const;
628
632 number
633 linfty_norm() const;
634
638 number
639 frobenius_norm() const;
640
645 number
646 trace() const;
647
653 void
654 invert();
655
664 void
665 solve(Vector<number> &v, const bool transposed = false) const;
666
671 void
672 solve(LAPACKFullMatrix<number> &B, const bool transposed = false) const;
673
692 void
693 compute_eigenvalues(const bool right_eigenvectors = false,
694 const bool left_eigenvectors = false);
695
715 void
716 compute_eigenvalues_symmetric(const number lower_bound,
717 const number upper_bound,
718 const number abs_accuracy,
721
748 void
751 const number lower_bound,
752 const number upper_bound,
753 const number abs_accuracy,
755 std::vector<Vector<number>> &eigenvectors,
756 const types::blas_int itype = 1);
757
773 void
776 std::vector<Vector<number>> &eigenvectors,
777 const types::blas_int itype = 1);
778
798 void
799 compute_svd();
800
820 void
821 compute_inverse_svd(const double threshold = 0.);
822
827 void
828 compute_inverse_svd_with_kernel(const unsigned int kernel_size);
829
836 std::complex<typename numbers::NumberTraits<number>::real_type>
837 eigenvalue(const size_type i) const;
838
850
857 get_left_eigenvectors() const;
858
863 number
864 singular_value(const size_type i) const;
865
870 inline const LAPACKFullMatrix<number> &
871 get_svd_u() const;
872
877 inline const LAPACKFullMatrix<number> &
878 get_svd_vt() const;
879
910 void
911 print_formatted(std::ostream &out,
912 const unsigned int precision = 3,
913 const bool scientific = true,
914 const unsigned int width = 0,
915 const char *zero_string = " ",
916 const double denominator = 1.,
917 const double threshold = 0.,
918 const char *separator = " ") const;
919
920private:
924 number
925 norm(const char type) const;
926
932
938
942 mutable std::vector<number> work;
943
947 mutable std::vector<types::blas_int> iwork;
948
955 std::vector<types::blas_int> ipiv;
956
960 std::vector<number> inv_work;
961
966 std::vector<typename numbers::NumberTraits<number>::real_type> wr;
967
972 std::vector<number> wi;
973
977 std::vector<number> vl;
978
982 std::vector<number> vr;
983
988 std::unique_ptr<LAPACKFullMatrix<number>> svd_u;
989
994 std::unique_ptr<LAPACKFullMatrix<number>> svd_vt;
995
1000};
1001
1002
1003
1009template <typename number>
1030
1031/*---------------------- Inline functions -----------------------------------*/
1032
1033template <typename number>
1034inline void
1036 const size_type j,
1037 const number value)
1038{
1039 (*this)(i, j) = value;
1040}
1041
1042
1043
1044template <typename number>
1047{
1048 return static_cast<size_type>(this->n_rows());
1049}
1050
1051
1052
1053template <typename number>
1056{
1057 return static_cast<size_type>(this->n_cols());
1058}
1059
1060
1061
1062template <typename number>
1063template <typename MatrixType>
1064inline void
1066{
1067 this->reinit(M.m(), M.n());
1068
1069 // loop over the elements of the argument matrix row by row, as suggested
1070 // in the documentation of the sparse matrix iterator class, and
1071 // copy them into the current object
1072 for (size_type row = 0; 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 this->el(row, entry->column()) = entry->value();
1079 }
1080
1081 state = LAPACKSupport::matrix;
1082}
1083
1084
1085
1086template <typename number>
1087template <typename MatrixType>
1088inline void
1090 const size_type dst_offset_i,
1091 const size_type dst_offset_j,
1092 const size_type src_offset_i,
1093 const size_type src_offset_j,
1094 const number factor,
1095 const bool transpose)
1096{
1097 // loop over the elements of the argument matrix row by row, as suggested
1098 // in the documentation of the sparse matrix iterator class
1099 for (size_type row = src_offset_i; row < M.m(); ++row)
1100 {
1101 const typename MatrixType::const_iterator end_row = M.end(row);
1102 for (typename MatrixType::const_iterator entry = M.begin(row);
1103 entry != end_row;
1104 ++entry)
1105 {
1106 const size_type i = transpose ? entry->column() : row;
1107 const size_type j = transpose ? row : entry->column();
1108
1109 const size_type dst_i = dst_offset_i + i - src_offset_i;
1110 const size_type dst_j = dst_offset_j + j - src_offset_j;
1111 if (dst_i < this->n_rows() && dst_j < this->n_cols())
1112 (*this)(dst_i, dst_j) = factor * entry->value();
1113 }
1114 }
1115
1116 state = LAPACKSupport::matrix;
1117}
1118
1119
1120
1121template <typename number>
1122template <typename number2>
1123void
1125 const Vector<number2> &,
1126 const bool) const
1127{
1128 Assert(false,
1129 ExcMessage("LAPACKFullMatrix<number>::vmult must be called with a "
1130 "matching Vector<double> vector type."));
1131}
1132
1133
1134
1135template <typename number>
1136template <typename number2>
1137void
1139 const Vector<number2> &) const
1140{
1141 Assert(false,
1142 ExcMessage("LAPACKFullMatrix<number>::vmult_add must be called with a "
1143 "matching Vector<double> vector type."));
1144}
1145
1146
1147
1148template <typename number>
1149template <typename number2>
1150void
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
1162template <typename number>
1163template <typename number2>
1164void
1166 const Vector<number2> &) const
1167{
1168 Assert(false,
1169 ExcMessage("LAPACKFullMatrix<number>::Tvmult_add must be called "
1170 "with a matching Vector<double> vector type."));
1171}
1172
1173
1174
1175namespace internal
1176{
1177 namespace LAPACKFullMatrixImplementation
1178 {
1179 template <typename RealNumber>
1180 std::complex<RealNumber>
1181 pack_complex(const RealNumber &real_part, const RealNumber &imaginary_part)
1182 {
1183 return std::complex<RealNumber>(real_part, imaginary_part);
1184 }
1185
1186 // The eigenvalues in LAPACKFullMatrix with complex-valued matrices are
1187 // contained in the 'wi' array, ignoring the 'wr' array.
1188 template <typename Number>
1189 std::complex<Number>
1190 pack_complex(const Number &, const std::complex<Number> &complex_number)
1191 {
1192 return complex_number;
1193 }
1194 } // namespace LAPACKFullMatrixImplementation
1195} // namespace internal
1196
1197
1198
1199template <typename number>
1200inline std::complex<typename numbers::NumberTraits<number>::real_type>
1202{
1204 Assert(wr.size() == this->n_rows(), ExcInternalError());
1205 Assert(wi.size() == this->n_rows(), ExcInternalError());
1206 AssertIndexRange(i, this->n_rows());
1207
1209}
1210
1211
1212
1213template <typename number>
1214inline number
1216{
1219 AssertIndexRange(i, wr.size());
1220
1221 return wr[i];
1222}
1223
1224
1225
1226template <typename number>
1227inline const LAPACKFullMatrix<number> &
1229{
1232
1233 return *svd_u;
1234}
1235
1236
1237
1238template <typename number>
1239inline const LAPACKFullMatrix<number> &
1241{
1244
1245 return *svd_vt;
1246}
1247
1248
1249
1251
1252#endif
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 > &)
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
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 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 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
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
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
const TableIndices< N > & size() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcInvalidState()
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcState(State 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.
std::complex< RealNumber > pack_complex(const RealNumber &real_part, const RealNumber &imaginary_part)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
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)