36 namespace LAPACKFullMatrixImplementation
59 static_assert(std::is_same_v<T, double> || std::is_same_v<T, float>,
60 "Only implemented for double and float");
61 Assert(
matrix.size() ==
static_cast<std::size_t
>(n_rows * n_rows),
68 Assert(
static_cast<std::size_t
>(n_rows * n_rows) <=
72 Assert(
static_cast<std::size_t
>(n_rows * n_rows) <=
76 static_cast<std::size_t
>(2 * n_rows) <=
real_work.size(),
114 std::is_same_v<T, double> || std::is_same_v<T, float>,
115 "Only implemented for std::complex<double> and std::complex<float>");
116 Assert(
matrix.size() ==
static_cast<std::size_t
>(n_rows * n_rows),
121 Assert(
static_cast<std::size_t
>(n_rows * n_rows) <=
125 Assert(
static_cast<std::size_t
>(n_rows * n_rows) <=
152 template <
typename T>
169 Assert(
static_cast<std::size_t
>(n_rows * n_cols) ==
matrix.size(),
196 template <
typename T>
205 std::vector<std::complex<T>> &work,
213 Assert(
static_cast<std::size_t
>(n_rows * n_cols) ==
matrix.size(),
245template <
typename number>
254template <
typename number>
263template <
typename number>
272template <
typename number>
278 property =
M.property;
284template <
typename number>
294template <
typename number>
303 (*
this)(i,
j) = copy(i,
j);
308template <
typename number>
311 const std::array<number, 3> &
csr,
324 const number t =
A(i,
j);
333 const number t =
A(
j, i);
342template <
typename number>
362 (*this)(i,
j) = copy(
ii,
jj);
369template <
typename number>
379template <
typename number>
380template <
typename number2>
386 for (
size_type i = 0; i < this->m(); ++i)
388 (*
this)(i,
j) =
M(i,
j);
397template <
typename number>
398template <
typename number2>
404 for (
size_type i = 0; i < this->m(); ++i)
406 (*
this)(i,
j) =
M.el(i,
j);
415template <
typename number>
422 if (this->n_elements() != 0)
423 this->reset_values();
431template <
typename number>
440 const char type =
'G';
441 const number
cfrom = 1.;
448 number *
values = this->values.data();
460template <
typename number>
471 const char type =
'G';
472 const number
cto = 1.;
479 number *
values = this->values.data();
491template <
typename number>
509 number *
values = this->values.data();
519 template <
typename number>
548 const std::array<number, 3>
csr =
554 const number t =
A(i,
k);
556 z(i) = -
csr[1] * t +
csr[0] * z(i);
591 const std::array<number, 3>
csr =
597 const number t =
A(i,
k);
599 z(i) = -
csr[1] * t +
csr[0] * z(i);
606 template <
typename number>
609 const std::complex<number> ,
610 const Vector<std::complex<number>> & )
618template <
typename number>
646 (*
this)(i,
j) = (*
this)(
j, i);
658template <
typename number>
666 const number
alpha = 1.;
668 const number
null = 0.;
679 const char diag =
'N';
680 const char trans =
'N';
691 &uplo, &
trans, &diag, &
N, this->values.data(), &
lda, w.data(), &
incx);
719 std::lock_guard<std::mutex> lock(mutex);
736 for (size_type i = 0; i < wr.size(); ++i)
754 std::lock_guard<std::mutex> lock(mutex);
771 for (size_type i = 0; i < wr.size(); ++i)
794template <
typename number>
802 const number
alpha = 1.;
804 const number
null = 0.;
815 const char diag =
'N';
816 const char trans =
'T';
827 &uplo, &
trans, &diag, &
N, this->values.data(), &
lda, w.data(), &
incx);
856 std::lock_guard<std::mutex> lock(mutex);
874 for (
size_type i = 0; i < wr.size(); ++i)
892 std::lock_guard<std::mutex> lock(mutex);
910 for (size_type i = 0; i < wr.size(); ++i)
933template <
typename number>
943template <
typename number>
953template <
typename number>
968 const number
alpha = 1.;
988template <
typename number>
1002 const number
alpha = 1.;
1024template <
typename number>
1051 std::lock_guard<std::mutex> lock(mutex);
1053 work.resize(
kk *
nn);
1062 work[
j *
kk + i] =
V(i) *
B(i,
j);
1066 const number
alpha = 1.;
1075 this->values.data(),
1086template <
typename number>
1095#ifdef DEAL_II_LAPACK_WITH_MKL
1096 const number
one = 1.;
1107template <
typename number>
1124template <
typename number>
1139 const number
alpha = 1.;
1149 this->values.data(),
1170 this->values.data(),
1182template <
typename number>
1196 const number
alpha = 1.;
1218template <
typename number>
1233 const number
alpha = 1.;
1243 this->values.data(),
1264 this->values.data(),
1276template <
typename number>
1290 const number
alpha = 1.;
1312template <
typename number>
1327 const number
alpha = 1.;
1336 this->values.data(),
1347template <
typename number>
1361 const number
alpha = 1.;
1383template <
typename number>
1392 number *
const values = this->values.data();
1407template <
typename number>
1416template <
typename number>
1420 const char type(
'O');
1426template <
typename number>
1430 const char type(
'I');
1436template <
typename number>
1440 const char type(
'F');
1446template <
typename number>
1450 std::lock_guard<std::mutex> lock(mutex);
1454 ExcMessage(
"norms can be called in matrix state only."));
1458 const number *
const values = this->values.data();
1463 (type ==
'I' || type ==
'O') ? std::max<types::blas_int>(1,
N) : 0;
1471 (type ==
'I') ? std::max<types::blas_int>(1,
M) : 0;
1479template <
typename number>
1485 ExcMessage(
"Trace can be called in matrix state only."));
1489 for (
size_type i = 0; i < this->m(); ++i)
1490 tr += (*
this)(i, i);
1497template <
typename number>
1510 number *
const values = this->values.data();
1524template <
typename number>
1528 std::lock_guard<std::mutex> lock(mutex);
1533 const number *
values = this->values.data();
1557template <
typename number>
1561 std::lock_guard<std::mutex> lock(mutex);
1567 const number *
const values = this->values.data();
1573 const char norm =
'1';
1574 const char diag =
'N';
1595template <
typename number>
1605 std::fill(wr.begin(), wr.end(), 0.);
1606 ipiv.resize(8 *
mm);
1608 svd_u = std::make_unique<LAPACKFullMatrix<number>>(
mm,
mm);
1609 svd_vt = std::make_unique<LAPACKFullMatrix<number>>(
nn,
nn);
1617 std::vector<typename numbers::NumberTraits<number>::real_type>
real_work;
1621 std::size_t min =
std::min(this->m(), this->n());
1622 std::size_t max =
std::max(this->m(), this->n());
1624 std::max(5 * min * min + 5 * min, 2 * max * min + 2 * min * min + min));
1672template <
typename number>
1683 for (
size_type i = 0; i < wr.size(); ++i)
1686 wr[i] =
one / wr[i];
1695template <
typename number>
1698 const unsigned int kernel_size)
1706 const unsigned int n_wr = wr.size();
1708 wr[i] =
one / wr[i];
1716template <
typename number>
1725 number *
const values = this->values.data();
1731 compute_lu_factorization();
1734 inv_work.resize(
mm);
1740 compute_cholesky_factorization();
1748 this->el(i,
j) = this->el(
j, i);
1759template <
typename number>
1767 const number *
const values = this->values.data();
1794 "The matrix has to be either factorized or triangular."));
1802template <
typename number>
1813 const number *
const values = this->values.data();
1862 "The matrix has to be either factorized or triangular."));
1870template <
typename number>
1894 for (
size_type i = 0; i < this->m(); ++i)
1896 (ipiv[i] ==
types::blas_int(i + 1) ? this->el(i, i) : -this->el(i, i));
1902template <
typename number>
1917 const char jobvr = (right) ?
V :
N;
1918 const char jobvl = (left) ?
V :
N;
1932 std::vector<typename numbers::NumberTraits<number>::real_type>
real_work;
1979 std::to_string(-
info) +
1981 " parameter had an illegal value."));
1988 "Lapack error in geev: the QR algorithm failed to compute "
1989 "all the eigenvalues, and no eigenvectors have been computed."));
2007 template <
typename RealNumber>
2010 const std::vector<RealNumber> &vr,
2011 const std::complex<RealNumber> &eigenvalue,
2013 unsigned int &index)
2015 const std::size_t n =
result.n();
2016 if (eigenvalue.imag() != 0.)
2018 for (std::size_t
j = 0;
j < n; ++
j)
2020 result(
j, index).real(vr[index * n +
j]);
2021 result(
j, index + 1).real(vr[index * n +
j]);
2022 result(
j, index).imag(vr[(index + 1) * n +
j]);
2023 result(
j, index + 1).imag(-vr[(index + 1) * n +
j]);
2032 for (
unsigned int j = 0;
j < n; ++
j)
2033 result(
j, index).real(vr[index * n +
j]);
2042 template <
typename ComplexNumber>
2045 const std::vector<ComplexNumber> &vr,
2048 unsigned int &index)
2050 const std::size_t n =
result.n();
2051 for (
unsigned int j = 0;
j < n; ++
j)
2061template <
typename number>
2067 ExcMessage(
"Right eigenvectors are not available! Did you "
2068 "set the associated flag in compute_eigenvalues()?"));
2073 for (
unsigned int i = 0; i < n();)
2081template <
typename number>
2087 ExcMessage(
"Left eigenvectors are not available! Did you "
2088 "set the associated flag in compute_eigenvalues()?"));
2093 for (
unsigned int i = 0; i < n();)
2101template <
typename number>
2104 const number lower_bound,
2105 const number upper_bound,
2117 number *
const values_A = this->values.data();
2121 const char *
const jobz(&
V);
2122 const char *
const uplo(&
U);
2123 const char *
const range(&
V);
2125 std::vector<types::blas_int> iwork(
static_cast<size_type>(5 *
nn));
2196 std::to_string(-
info) +
2198 " parameter had an illegal value."));
2204 "Lapack error in syevx: " + std::to_string(
info) +
2205 " eigenvectors failed to converge."
2206 " (You may need to scale the abs_accuracy according"
2207 " to your matrix norm.)"));
2212 ExcMessage(
"Lapack error in syevx: unknown error."));
2233template <
typename number>
2237 const number lower_bound,
2238 const number upper_bound,
2253 number *
const values_A = this->values.data();
2258 const char *
const jobz(&
V);
2259 const char *
const uplo(&
U);
2260 const char *
const range(&
V);
2341 std::to_string(-
info) +
2343 " parameter had an illegal value."));
2350 "Lapack error in sygvx: ssyevx/dsyevx failed to converge, and " +
2351 std::to_string(
info) +
2352 " eigenvectors failed to converge."
2353 " (You may need to scale the abs_accuracy"
2354 " according to the norms of matrices A and B.)"));
2360 "Lapack error in sygvx: the leading minor of order " +
2361 std::to_string(
info -
nn) +
2362 " of matrix B is not positive-definite."
2363 " The factorization of B could not be completed and"
2364 " no eigenvalues or eigenvectors were computed."));
2369 ExcMessage(
"Lapack error in sygvx: unknown error."));
2391template <
typename number>
2404 ExcMessage(
"eigenvectors.size() > matrix.n()"));
2410 number *
const values_A = this->values.data();
2416 const char *
const uplo = (&
U);
2471 std::to_string(-
info) +
2473 " parameter had an illegal value."));
2480 "Lapack error in sygv: ssyev/dsyev failed to converge, and " +
2481 std::to_string(
info) +
2482 " off-diagonal elements of an intermediate "
2483 " tridiagonal did not converge to zero."
2484 " (You may need to scale the abs_accuracy"
2485 " according to the norms of matrices A and B.)"));
2491 "Lapack error in sygv: the leading minor of order " +
2492 std::to_string(
info -
nn) +
2493 " of matrix B is not positive-definite."
2494 " The factorization of B could not be completed and"
2495 " no eigenvalues or eigenvectors were computed."));
2500 ExcMessage(
"Lapack error in sygv: unknown error."));
2517template <
typename number>
2520 const unsigned int precision,
2521 const bool scientific,
2522 const unsigned int width_,
2525 const double threshold,
2526 const char *separator)
const
2528 unsigned int width =
width_;
2538 std::ios::fmtflags
old_flags = out.flags();
2543 out.setf(std::ios::scientific, std::ios::floatfield);
2545 width = precision + 7;
2549 out.setf(std::ios::fixed, std::ios::floatfield);
2551 width = precision + 2;
2554 for (
size_type i = 0; i < this->m(); ++i)
2562 out << std::setw(width) << (*this)(i,
j) << separator;
2563 else if (
std::abs(this->el(i,
j)) > threshold)
2564 out << std::setw(width) << this->el(i,
j) *
denominator << separator;
2566 out << std::setw(width) <<
zero_string << separator;
2578template <
typename number>
2588template <
typename number>
2597template <
typename number>
2607template <
typename number>
2613 matrix->solve(dst,
false);
2617template <
typename number>
2623 matrix->solve(dst,
true);
2627template <
typename number>
2635 matrix->solve(*aux,
false);
2640template <
typename number>
2648 matrix->solve(*aux,
true);
2654#include "lapack_full_matrix.inst"
void omatcopy(char, char, ::types::blas_int, ::types::blas_int, const number1, const number2 *, ::types::blas_int, number3 *, ::types::blas_int)
EnableObserverPointer & operator=(const EnableObserverPointer &)
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 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
LAPACKSupport::State get_state() const
void compute_eigenvalues_symmetric(const number lower_bound, const number upper_bound, const number abs_accuracy, Vector< number > &eigenvalues, FullMatrix< number > &eigenvectors)
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
void reinit(const size_type size)
void compute_cholesky_factorization()
LAPACKFullMatrix< number > & operator=(const LAPACKFullMatrix< number > &)
void compute_lu_factorization()
FullMatrix< std::complex< typename numbers::NumberTraits< number >::real_type > > get_left_eigenvectors() const
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)
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)
std::make_unsigned_t< types::blas_int > size_type
number frobenius_norm() const
LAPACKFullMatrix(const size_type size=0)
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
number linfty_norm() const
void TmTmult(LAPACKFullMatrix< number > &C, const LAPACKFullMatrix< number > &B, const bool adding=false) const
void compute_inverse_svd_with_kernel(const unsigned int kernel_size)
void Tvmult_add(Vector< number2 > &w, const Vector< number2 > &v) const
void rank1_update(const number a, const Vector< number > &v)
void remove_row_and_column(const size_type row, const size_type col)
LAPACKFullMatrix< number > & operator/=(const number factor)
number determinant() const
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 vmult(Vector< number > &, const Vector< number > &) const
void initialize(const LAPACKFullMatrix< number > &)
void Tvmult(Vector< number > &, const Vector< number > &) const
numbers::NumberTraits< Number >::real_type norm() const
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
void reinit(const size_type size1, const size_type size2, const bool omit_default_initialization=false)
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcZero()
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcErrorCode(std::string arg1, types::blas_int arg2)
static ::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
static ::ExceptionBase & ExcProperty(Property arg1)
#define Assert(cond, exc)
static ::ExceptionBase & ExcSingular()
#define AssertIsFinite(number)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcNotQuadratic()
static ::ExceptionBase & ExcInvalidState()
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcState(State arg1)
#define AssertThrow(cond, exc)
void getrs(const char *, const ::types::blas_int *, const ::types::blas_int *, const number1 *, const ::types::blas_int *, const ::types::blas_int *, number2 *, const ::types::blas_int *, ::types::blas_int *)
void syrk(const char *, const char *, const ::types::blas_int *, const ::types::blas_int *, const number1 *, const number2 *, const ::types::blas_int *, const number3 *, number4 *, const ::types::blas_int *)
void geev(const char *, const char *, const ::types::blas_int *, number1 *, const ::types::blas_int *, number2 *, number3 *, number4 *, const ::types::blas_int *, number5 *, const ::types::blas_int *, number6 *, const ::types::blas_int *, ::types::blas_int *)
void gemm(const char *, const char *, const ::types::blas_int *, const ::types::blas_int *, const ::types::blas_int *, const number1 *, const number2 *, const ::types::blas_int *, const number3 *, const ::types::blas_int *, const number4 *, number5 *, const ::types::blas_int *)
void pocon(const char *, const ::types::blas_int *, const number1 *, const ::types::blas_int *, const number2 *, number3 *, number4 *, ::types::blas_int *, ::types::blas_int *)
void lascl(const char *, const ::types::blas_int *, const ::types::blas_int *, const number1 *, const number2 *, const ::types::blas_int *, const ::types::blas_int *, number3 *, const ::types::blas_int *, ::types::blas_int *)
void syr(const char *, const ::types::blas_int *, const number1 *, const number2 *, const ::types::blas_int *, number3 *, const ::types::blas_int *)
void trtrs(const char *, const char *, const char *, const ::types::blas_int *, const ::types::blas_int *, const number1 *, const ::types::blas_int *, number2 *, const ::types::blas_int *, ::types::blas_int *)
void syevx(const char *, const char *, const char *, const ::types::blas_int *, number1 *, const ::types::blas_int *, const number2 *, const number3 *, const ::types::blas_int *, const ::types::blas_int *, const number4 *, ::types::blas_int *, number5 *, number6 *, const ::types::blas_int *, number7 *, const ::types::blas_int *, ::types::blas_int *, ::types::blas_int *, ::types::blas_int *)
void trcon(const char *, const char *, const char *, const ::types::blas_int *, const number1 *, const ::types::blas_int *, number2 *, number3 *, ::types::blas_int *, ::types::blas_int *)
void gesdd(const char *, const ::types::blas_int *, const ::types::blas_int *, number1 *, const ::types::blas_int *, number2 *, number3 *, const ::types::blas_int *, number4 *, const ::types::blas_int *, number5 *, const ::types::blas_int *, ::types::blas_int *, ::types::blas_int *)
void axpy(const ::types::blas_int *, const number1 *, const number2 *, const ::types::blas_int *, number3 *, const ::types::blas_int *)
void trmv(const char *, const char *, const char *, const ::types::blas_int *, const number1 *, const ::types::blas_int *, number2 *, const ::types::blas_int *)
void sygvx(const ::types::blas_int *, const char *, const char *, const char *, const ::types::blas_int *, number1 *, const ::types::blas_int *, number2 *, const ::types::blas_int *, const number3 *, const number4 *, const ::types::blas_int *, const ::types::blas_int *, const number5 *, ::types::blas_int *, number6 *, number7 *, const ::types::blas_int *, number8 *, const ::types::blas_int *, ::types::blas_int *, ::types::blas_int *, ::types::blas_int *)
void gemv(const char *, const ::types::blas_int *, const ::types::blas_int *, const number1 *, const number2 *, const ::types::blas_int *, const number3 *, const ::types::blas_int *, const number4 *, number5 *, const ::types::blas_int *)
number1 lange(const char *, const ::types::blas_int *, const ::types::blas_int *, const number1 *, const ::types::blas_int *, number2 *)
void potrs(const char *, const ::types::blas_int *, const ::types::blas_int *, const number1 *, const ::types::blas_int *, number2 *, const ::types::blas_int *, ::types::blas_int *)
void sygv(const ::types::blas_int *, const char *, const char *, const ::types::blas_int *, number1 *, const ::types::blas_int *, number2 *, const ::types::blas_int *, number3 *, number4 *, const ::types::blas_int *, ::types::blas_int *)
void potri(const char *, const ::types::blas_int *, number1 *, const ::types::blas_int *, ::types::blas_int *)
number1 lansy(const char *, const char *, const ::types::blas_int *, const number1 *, const ::types::blas_int *, number2 *)
void potrf(const char *, const ::types::blas_int *, number1 *, const ::types::blas_int *, ::types::blas_int *)
void getrf(const ::types::blas_int *, const ::types::blas_int *, number1 *, const ::types::blas_int *, ::types::blas_int *, ::types::blas_int *)
void getri(const ::types::blas_int *, number1 *, const ::types::blas_int *, const ::types::blas_int *, number2 *, const ::types::blas_int *, ::types::blas_int *)
@ cholesky
Contents is a Cholesky decomposition.
@ lu
Contents is an LU decomposition.
@ matrix
Contents is actually a matrix.
@ unusable
Contents is something useless.
@ inverse_matrix
Contents is the inverse of a matrix.
@ svd
Matrix contains singular value decomposition,.
@ inverse_svd
Matrix is the inverse of a singular value decomposition.
@ eigenvalues
Eigenvalue vector is filled.
@ symmetric
Matrix is symmetric.
@ upper_triangular
Matrix is upper triangular.
@ lower_triangular
Matrix is lower triangular.
@ general
No special properties.
static const types::blas_int one
std::array< NumberType, 3 > givens_rotation(const NumberType &x, const NumberType &y)
std::array< NumberType, 3 > hyperbolic_rotation(const NumberType &x, const NumberType &y)
void gesdd_helper(const char job, const types::blas_int n_rows, const types::blas_int n_cols, AlignedVector< T > &matrix, std::vector< T > &singular_values, AlignedVector< T > &left_vectors, AlignedVector< T > &right_vectors, std::vector< T > &real_work, std::vector< T > &, std::vector< types::blas_int > &integer_work, const types::blas_int work_flag, types::blas_int &info)
void geev_helper(const char vl, const char vr, AlignedVector< T > &matrix, const types::blas_int n_rows, std::vector< T > &real_part_eigenvalues, std::vector< T > &imag_part_eigenvalues, std::vector< T > &left_eigenvectors, std::vector< T > &right_eigenvectors, std::vector< T > &real_work, std::vector< T > &, const types::blas_int work_flag, types::blas_int &info)
bool is_nan(const double x)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
static bool equal(const T *p1, const T *p2)
static constexpr const number & conjugate(const number &x)
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)