Reference documentation for deal.II version 9.5.0
\(\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// 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
37template <typename number>
38class Vector;
39template <typename number>
40class BlockVector;
41template <typename number>
42class FullMatrix;
43template <typename number>
44class SparseMatrix;
45#endif
46
58template <typename number>
59class LAPACKFullMatrix : public TransposeTable<number>
60{
61public:
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
281 m() const;
282
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
431 const bool adding = false) const;
432
437 void
440 const bool adding = false) const;
441
456 void
459 const bool adding = false) const;
460
465 void
468 const bool adding = false) const;
469
486 void
489 const Vector<number> & V,
490 const bool adding = false) const;
491
506 void
509 const bool adding = false) const;
510
515 void
518 const bool adding = false) const;
519
535 void
538 const bool adding = false) const;
539
544 void
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
717 compute_eigenvalues_symmetric(const number lower_bound,
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
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
918private:
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
1007template <typename number>
1009{
1010public:
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
1021 void
1023
1024private:
1027};
1028
1029/*---------------------- Inline functions -----------------------------------*/
1030
1031template <typename number>
1032inline void
1034 const size_type j,
1035 const number value)
1036{
1037 (*this)(i, j) = value;
1038}
1039
1040
1041
1042template <typename number>
1045{
1046 return static_cast<size_type>(this->n_rows());
1047}
1048
1049
1050
1051template <typename number>
1054{
1055 return static_cast<size_type>(this->n_cols());
1056}
1057
1058
1059
1060template <typename number>
1061template <typename MatrixType>
1062inline 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
1084template <typename number>
1085template <typename MatrixType>
1086inline 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
1119template <typename number>
1120template <typename number2>
1121void
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
1133template <typename number>
1134template <typename number2>
1135void
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
1146template <typename number>
1147template <typename number2>
1148void
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
1160template <typename number>
1161template <typename number2>
1162void
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
1173namespace 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
1197template <typename number>
1198inline 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
1211template <typename number>
1212inline number
1214{
1217 AssertIndexRange(i, wr.size());
1218
1219 return wr[i];
1220}
1221
1222
1223
1224template <typename number>
1225inline const LAPACKFullMatrix<number> &
1227{
1230
1231 return *svd_u;
1232}
1233
1234
1235
1236template <typename number>
1237inline const LAPACKFullMatrix<number> &
1239{
1242
1243 return *svd_vt;
1244}
1245
1246
1247
1249
1250#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)
std::make_unsigned< types::blas_int >::type size_type
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
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
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 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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
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)