*/
explicit LAPACKFullMatrix(const size_type size = 0);
-
/**
* Constructor. Initialize the matrix as a rectangular matrix $\rm{rows}
* \times \rm{cols}$.
*/
LAPACKFullMatrix(const size_type rows, const size_type cols);
-
/**
* Copy constructor. This constructor does a deep copy of the matrix.
* Therefore, it poses a possible efficiency problem, if for example,
* @note It is assumed that @p A, @p B and @p V have compatible sizes and that
* @p C already has the right size.
*
- * @note This function is not provided by LAPACK. The function first forms $\rm{diag}(\mathbf V) \cdot \mathbf B$ product and
- * then uses Xgemm function.
+ * @note This function is not provided by LAPACK. The function first forms
+ * $\rm{diag}(\mathbf V) \cdot \mathbf B$ product and then uses the Xgemm
+ * function.
*/
void
Tmmult(LAPACKFullMatrix<number> & C,
compute_inverse_svd_with_kernel(const unsigned int kernel_size);
/**
- * Retrieve eigenvalue after compute_eigenvalues() was called.
+ * Retrieve eigenvalue after compute_eigenvalues() was called. The return
+ * type is always a complex number: In case the underlying matrix is
+ * real-valued, the type `std::complex<number>` is returned, whereas
+ * `number` is returned for complex numbers.
*/
- std::complex<number>
+ std::complex<typename numbers::NumberTraits<number>::real_type>
eigenvalue(const size_type i) const;
/**
- * After a call to compute_eigenvalues(), this function returns all (right)
- * eigenvectors as returned by LAPACK. This means that eigenvectors are
- * contained in column-major ordering in a matrix associated to the given
- * flat vector. Note that for real-valued matrices, there might appear
- * complex eigenvalues with complex-conjugate values and eigenvectors. For
- * those cases, LAPACK places a column of n() entries with the real part and
- * the next column (corresponding to the complex conjugate eigenvalue) for
- * the imaginary part of the eigenvector.
+ * After a call to compute_eigenvalues(), this function returns the $n\times
+ * n$ matrix of (right) eigenvectors in a decomposition of the form $A V = V
+ * \Lambda$. Note that this function constructs the associated matrix on the
+ * fly, since LAPACK packs complex-conjugate eigenvalue/eigenvector pairs of
+ * real-valued matrices into a real-valued return matrix. This call only
+ * succeeds in case the respective flag `right_eigenvectors` in
+ * compute_eigenvalues() has been set to true.
*/
- std::vector<number>
+ FullMatrix<std::complex<typename numbers::NumberTraits<number>::real_type>>
get_right_eigenvectors() const;
/**
* Return the matrix of left eigenvectors after a call to
- * compute_eigenvalues(), following the same convention as
+ * compute_eigenvalues(), following the same principles as
* get_right_eigenvectors().
*/
- std::vector<number>
+ FullMatrix<std::complex<typename numbers::NumberTraits<number>::real_type>>
get_left_eigenvectors() const;
/**
}
+
template <typename number>
inline typename LAPACKFullMatrix<number>::size_type
LAPACKFullMatrix<number>::m() const
return static_cast<size_type>(this->n_rows());
}
+
+
template <typename number>
inline typename LAPACKFullMatrix<number>::size_type
LAPACKFullMatrix<number>::n() const
return static_cast<size_type>(this->n_cols());
}
+
+
template <typename number>
template <typename MatrixType>
inline void
}
+
template <typename number>
template <typename number2>
void
}
+
template <typename number>
template <typename number2>
void
}
+
template <typename number>
template <typename number2>
void
}
+
template <typename number>
template <typename number2>
void
const Vector<number2> &) const
{
Assert(false,
- ExcMessage(
- "LAPACKFullMatrix<number>::Tvmult_add must be called with a "
- "matching Vector<double> vector type."));
+ ExcMessage("LAPACKFullMatrix<number>::Tvmult_add must be called "
+ "with a matching Vector<double> vector type."));
}
-template <typename number>
-inline std::complex<number>
-LAPACKFullMatrix<number>::eigenvalue(const size_type i) const
-{
- Assert(state & LAPACKSupport::eigenvalues, ExcInvalidState());
- Assert(wr.size() == this->n_rows(), ExcInternalError());
- Assert(wi.size() == this->n_rows(), ExcInternalError());
- AssertIndexRange(i, this->n_rows());
-
- if (numbers::NumberTraits<number>::is_complex)
- return std::complex<number>(wi[i]);
- else
- return std::complex<number>(wr[i], wi[i]);
-}
-
-template <typename number>
-inline std::vector<number>
-LAPACKFullMatrix<number>::get_right_eigenvectors() const
+namespace internal
{
- Assert(state & LAPACKSupport::eigenvalues, ExcInvalidState());
- Assert(vr.size() == this->n_rows() * this->n_cols(),
- ExcMessage("Right eigenvectors are not available! Did you "
- "set the associated flag in compute_eigenvalues()"));
+ namespace LAPACKFullMatrixImplementation
+ {
+ template <typename RealNumber>
+ std::complex<RealNumber>
+ pack_complex(const RealNumber &real_part, const RealNumber &imaginary_part)
+ {
+ return std::complex<RealNumber>(real_part, imaginary_part);
+ }
+
+ // The eigenvalues in LAPACKFullMatrix with complex-valued matrices are
+ // contained in the 'wi' array, ignoring the 'wr' array.
+ template <typename Number>
+ std::complex<Number>
+ pack_complex(const Number &, const std::complex<Number> &complex_number)
+ {
+ return complex_number;
+ }
+ } // namespace LAPACKFullMatrixImplementation
+} // namespace internal
- return vr;
-}
template <typename number>
-inline std::vector<number>
-LAPACKFullMatrix<number>::get_left_eigenvectors() const
+inline std::complex<typename numbers::NumberTraits<number>::real_type>
+LAPACKFullMatrix<number>::eigenvalue(const size_type i) const
{
Assert(state & LAPACKSupport::eigenvalues, ExcInvalidState());
- Assert(vl.size() == this->n_rows() * this->n_cols(),
- ExcMessage("Left eigenvectors are not available! Did you "
- "set the associated flag in compute_eigenvalues()"));
+ Assert(wr.size() == this->n_rows(), ExcInternalError());
+ Assert(wi.size() == this->n_rows(), ExcInternalError());
+ AssertIndexRange(i, this->n_rows());
- return vl;
+ return internal::LAPACKFullMatrixImplementation::pack_complex(wr[i], wi[i]);
}
+
template <typename number>
inline number
LAPACKFullMatrix<number>::singular_value(const size_type i) const
}
+
template <typename number>
inline const LAPACKFullMatrix<number> &
LAPACKFullMatrix<number>::get_svd_u() const
}
+
template <typename number>
inline const LAPACKFullMatrix<number> &
LAPACKFullMatrix<number>::get_svd_vt() const
&info);
}
+
+
template <typename T>
void
geev_helper(const char vl,
Assert(static_cast<std::size_t>(n_rows * n_rows) <=
right_eigenvectors.size(),
ExcInternalError());
- Assert(std::max<std::size_t>(1, work_flag) <= real_work.size(),
- ExcInternalError());
Assert(work_flag == -1 ||
- std::max<long int>(1, 2 * n_rows) <= (work_flag),
+ std::max<std::size_t>(1, work_flag) <= real_work.size(),
+ ExcInternalError());
+ Assert(work_flag == -1 || std::max<long int>(1, 2 * n_rows) <= work_flag,
ExcInternalError());
geev(&vl,
} // namespace LAPACKFullMatrixImplementation
} // namespace internal
+
+
template <typename number>
LAPACKFullMatrix<number>::LAPACKFullMatrix(const size_type n)
: TransposeTable<number>(n, n)
{}
+
template <typename number>
LAPACKFullMatrix<number>::LAPACKFullMatrix(const size_type m, const size_type n)
: TransposeTable<number>(m, n)
{}
+
template <typename number>
LAPACKFullMatrix<number>::LAPACKFullMatrix(const LAPACKFullMatrix &M)
: TransposeTable<number>(M)
{}
+
template <typename number>
LAPACKFullMatrix<number> &
LAPACKFullMatrix<number>::operator=(const LAPACKFullMatrix<number> &M)
}
+
template <typename number>
template <typename number2>
LAPACKFullMatrix<number> &
}
+
template <typename number>
template <typename number2>
LAPACKFullMatrix<number> &
}
+
template <typename number>
LAPACKFullMatrix<number> &
LAPACKFullMatrix<number>::operator=(const number d)
}
+
template <typename number>
LAPACKFullMatrix<number> &
LAPACKFullMatrix<number>::operator*=(const number factor)
}
+
template <typename number>
LAPACKFullMatrix<number> &
LAPACKFullMatrix<number>::operator/=(const number factor)
}
+
template <typename number>
void
LAPACKFullMatrix<number>::Tvmult(Vector<number> & w,
}
+
template <typename number>
void
LAPACKFullMatrix<number>::vmult_add(Vector<number> & w,
}
+
template <typename number>
void
LAPACKFullMatrix<number>::Tvmult_add(Vector<number> & w,
}
+
template <typename number>
void
LAPACKFullMatrix<number>::mmult(LAPACKFullMatrix<number> & C,
}
+
template <typename number>
void
LAPACKFullMatrix<number>::mmult(FullMatrix<number> & C,
}
+
template <typename number>
void
LAPACKFullMatrix<number>::scale_rows(const Vector<number> &V)
}
+
template <typename number>
void
LAPACKFullMatrix<number>::Tmmult(FullMatrix<number> & C,
}
+
template <typename number>
void
LAPACKFullMatrix<number>::mTmult(LAPACKFullMatrix<number> & C,
}
+
template <typename number>
void
LAPACKFullMatrix<number>::TmTmult(LAPACKFullMatrix<number> & C,
}
+
template <typename number>
void
LAPACKFullMatrix<number>::TmTmult(FullMatrix<number> & C,
}
+
template <typename number>
void
LAPACKFullMatrix<number>::compute_lu_factorization()
}
+
template <typename number>
void
LAPACKFullMatrix<number>::compute_inverse_svd(const double threshold)
}
+
template <typename number>
void
LAPACKFullMatrix<number>::compute_eigenvalues(const bool right, const bool left)
// resize workspace array
work.resize(lwork);
+ real_work.resize(lwork);
// Finally compute the eigenvalues.
internal::LAPACKFullMatrixImplementation::geev_helper(jobvl,
}
+
+namespace
+{
+ // This function extracts complex eigenvectors from the underlying 'number'
+ // array 'vr' of the LAPACK eigenvalue routine. For real-valued matrices
+ // addressed by this function specialization, we might get complex
+ // eigenvalues, which come in complex-conjugate pairs. In LAPACK, a compact
+ // storage scheme is applied that stores the real and imaginary part of
+ // eigenvectors only once, putting the real parts in one column and the
+ // imaginary part in the next of a real-valued array. Here, we do the
+ // unpacking into the usual complex values.
+ template <typename RealNumber>
+ void
+ unpack_lapack_eigenvector_and_increment_index(
+ const std::vector<RealNumber> & vr,
+ const std::complex<RealNumber> & eigenvalue,
+ FullMatrix<std::complex<RealNumber>> &result,
+ unsigned int & index)
+ {
+ const std::size_t n = result.n();
+ if (eigenvalue.imag() != 0.)
+ {
+ for (std::size_t j = 0; j < n; ++j)
+ {
+ result(j, index).real(vr[index * n + j]);
+ result(j, index + 1).real(vr[index * n + j]);
+ result(j, index).imag(vr[(index + 1) * n + j]);
+ result(j, index + 1).imag(-vr[(index + 1) * n + j]);
+ }
+
+ // we filled two columns with the complex-conjugate pair, so increment
+ // returned index by 2
+ index += 2;
+ }
+ else
+ {
+ for (unsigned int j = 0; j < n; ++j)
+ result(j, index).real(vr[index * n + j]);
+
+ // real-valued case, we only filled one column
+ ++index;
+ }
+ }
+
+ // This specialization fills the eigenvectors for complex-valued matrices,
+ // in which case we simply read off the entry in the 'vr' array.
+ template <typename ComplexNumber>
+ void
+ unpack_lapack_eigenvector_and_increment_index(
+ const std::vector<ComplexNumber> &vr,
+ const ComplexNumber &,
+ FullMatrix<ComplexNumber> &result,
+ unsigned int & index)
+ {
+ const std::size_t n = result.n();
+ for (unsigned int j = 0; j < n; ++j)
+ result(j, index) = vr[index * n + j];
+
+ // complex-valued case always only fills a single column
+ ++index;
+ }
+} // namespace
+
+
+
+template <typename number>
+FullMatrix<std::complex<typename numbers::NumberTraits<number>::real_type>>
+LAPACKFullMatrix<number>::get_right_eigenvectors() const
+{
+ Assert(state & LAPACKSupport::eigenvalues, ExcInvalidState());
+ Assert(vr.size() == this->n_rows() * this->n_cols(),
+ ExcMessage("Right eigenvectors are not available! Did you "
+ "set the associated flag in compute_eigenvalues()?"));
+
+ FullMatrix<std::complex<typename numbers::NumberTraits<number>::real_type>>
+ result(n(), n());
+
+ for (unsigned int i = 0; i < n();)
+ unpack_lapack_eigenvector_and_increment_index(vr, eigenvalue(i), result, i);
+
+ return result;
+}
+
+
+
+template <typename number>
+FullMatrix<std::complex<typename numbers::NumberTraits<number>::real_type>>
+LAPACKFullMatrix<number>::get_left_eigenvectors() const
+{
+ Assert(state & LAPACKSupport::eigenvalues, ExcInvalidState());
+ Assert(vl.size() == this->n_rows() * this->n_cols(),
+ ExcMessage("Left eigenvectors are not available! Did you "
+ "set the associated flag in compute_eigenvalues()?"));
+
+ FullMatrix<std::complex<typename numbers::NumberTraits<number>::real_type>>
+ result(n(), n());
+
+ for (unsigned int i = 0; i < n();)
+ unpack_lapack_eigenvector_and_increment_index(vl, eigenvalue(i), result, i);
+
+ return result;
+}
+
+
+
template <typename number>
void
LAPACKFullMatrix<number>::compute_eigenvalues_symmetric(
}
+
template <typename number>
void
LAPACKFullMatrix<number>::compute_generalized_eigenvalues_symmetric(
}
+
template <typename number>
void
LAPACKFullMatrix<number>::compute_generalized_eigenvalues_symmetric(
}
+
template <typename number>
void
LAPACKFullMatrix<number>::print_formatted(std::ostream & out,
deallog << lambda << " ";
}
deallog << std::endl;
- deallog << "Right eigenvectors ";
- for (const auto d : LA.get_right_eigenvectors())
- deallog << d << " ";
- deallog << std::endl;
- deallog << "Left eigenvectors ";
- for (const auto d : LA.get_left_eigenvectors())
- deallog << d << " ";
+ deallog << "Right eigenvectors" << std::endl;
+ LA.get_right_eigenvectors().print_formatted(deallog.get_file_stream());
+ deallog << "Left eigenvectors" << std::endl;
+ LA.get_left_eigenvectors().print_formatted(deallog.get_file_stream());
deallog << std::endl;
}
4.000e+00 -1.000e+00
-1.000e+00 4.000e+00
DEAL::Eigenvalues: (5.00000,0.00000) (3.00000,0.00000)
-DEAL::Right eigenvectors 0.707107 -0.707107 0.707107 0.707107
-DEAL::Left eigenvectors 0.707107 -0.707107 0.707107 0.707107
+DEAL::Right eigenvectors
+(7.071e-01,0.000e+00) (7.071e-01,0.000e+00)
+(-7.071e-01,0.000e+00) (7.071e-01,0.000e+00)
+DEAL::Left eigenvectors
+(7.071e-01,0.000e+00) (7.071e-01,0.000e+00)
+(-7.071e-01,0.000e+00) (7.071e-01,0.000e+00)
+DEAL::
DEAL::Checking matrix
1.000e+00 -1.000e+00
1.000e+00 1.000e+00
DEAL::Eigenvalues: (1.00000,1.00000) (1.00000,-1.00000)
-DEAL::Right eigenvectors 0.707107 0.00000 0.00000 -0.707107
-DEAL::Left eigenvectors -0.707107 0.00000 0.00000 0.707107
+DEAL::Right eigenvectors
+(7.071e-01,0.000e+00) (7.071e-01,0.000e+00)
+(0.000e+00,-7.071e-01) (0.000e+00,7.071e-01)
+DEAL::Left eigenvectors
+(-7.071e-01,0.000e+00) (-7.071e-01,0.000e+00)
+(0.000e+00,7.071e-01) (0.000e+00,-7.071e-01)
+DEAL::
projects / dealii.git / commitdiff