Line data    Source code

       1            1 : // Distributed under the MIT License.
       2              : // See LICENSE.txt for details.
       3              : 
       4              : /// \file
       5              : /// Defines function find_real_eigensystem.
       6              : 
       7              : #pragma once
       8              : 
       9              : #include <cstddef>
      10              : 
      11              : #include "DataStructures/Tensor/TypeAliases.hpp"
      12              : 
      13              : /// \cond
      14              : class DataVector;
      15              : namespace gsl {
      16              : template <typename T>
      17              : class not_null;
      18              : }  // namespace gsl
      19              : /// \endcond
      20              : 
      21              : /*!
      22              :  * \ingroup NumericalAlgorithmsGroup
      23              :  * \brief Find the eigensystem of a diagonalizable matrix whose eigenvalues are
      24              :  * all real, with a complete set of eigenvectors even when degeneracies are
      25              :  * present.
      26              :  *
      27              :  * We label the matrix $A_{m}{}^{n}$ such that the index $m$ labels rows and the
      28              :  * index $n$ labels columns.
      29              :  *
      30              :  * With this choice of indices, the $i$th left eigensystem is given by
      31              :  * \begin{equation}
      32              :  *   \sum_m L^{im} A_{m}{}^{n} = \lambda_{(i)} L^{in},
      33              :  * \end{equation}
      34              :  * where $\lambda_{(i)}$ is the $i$th eigenvalue.
      35              :  * Similarly, the $i$th right eigensystem is given by
      36              :  * \begin{equation}
      37              :  *   \sum_n R_{in} A_{m}{}^{n} = \lambda_{(i)} R_{im}.
      38              :  * \end{equation}
      39              :  *
      40              :  * We take the matrix and provide the eigensystem in the following types:
      41              :  * - %Matrix $A_{m}{}^{n}$: `tnsr::iJ<DataVector, MatrixSize> matrix`
      42              :  * - Eigenvalues $\lambda_{(i)}$: `tnsr::i<DataVector, MatrixSize> eigenvalues`
      43              :  * - Right eigenvectors $R_{in}$: `tnsr::ij<DataVector, MatrixSize>
      44              :  *   right_eigenvectors`
      45              :  * - Left eigenvectors $L^{im}$: `tnsr::IJ<DataVector, MatrixSize>
      46              :  *   left_eigenvectors`
      47              :  *
      48              :  * Note that the indices $i$, $m$, and $n$ are not tensorial. They simply label
      49              :  * the eigensystem, rows of the matrix, and columns of the matrix, respectively.
      50              :  * We make the conventions (raised or lowered) presented above to clearly
      51              :  * distinguish which indices (first or second) of the matrix get contracted
      52              :  * with which eigenvectors (left or right).
      53              :  *
      54              :  * Each grid point holds an independent matrix. Repeated (degenerate)
      55              :  * eigenvalues are allowed, as long as the matrix is diagonalizable. The
      56              :  * eigenvalues are in increasing order, with degenerate eigenvalues grouped
      57              :  * together.
      58              :  *
      59              :  * To compute such an eigensystem, we follow the procedure outlined below, each
      60              :  * step using a LAPACK routine.
      61              :  *
      62              :  * 1. Balance the matrix (`dgebal`), $A' = D^{-1} A D$, where $D$ is a diagonal
      63              :  *    matrix chosen so that the rows and columns of $A'$ have similar norms.
      64              :  *    This improves the accuracy of the next steps and does not change the
      65              :  *    eigenvalues.
      66              :  * 2. Compute the real Schur form of the balanced matrix (`dgees`),
      67              :  *    $A' = Z T Z^T$, where $Z$ is orthogonal and $T$ is upper triangular with
      68              :  *    the eigenvalues on its diagonal. The first column of $Z$ is a right
      69              :  *    eigenvector of $A'$ for the first eigenvalue on the diagonal of $T$. In
      70              :  *    general, the other columns are not eigenvectors. The exception is a
      71              :  *    repeated eigenvalue (i.e., when we have a degeneracy): if the first $k$
      72              :  *    eigenvalues on the diagonal of $T$ are equal, the first $k$ columns of $Z$
      73              :  *    are $k$ orthonormal right eigenvectors for that degenerate subspace.
      74              :  * 3. Group the eigenvalues into degenerate families, i.e., sets of equal
      75              :  *    eigenvalues. For each family, of size $k$, reorder the Schur form
      76              :  *    (`dtrsen`), $T \to Q^T T Q$ and $Z \to Z Q$ with $Q$ orthogonal, so that
      77              :  *    the family's eigenvalues come first on the diagonal of $T$, and take the
      78              :  *    first $k$ columns of $Z$ as the family's right eigenvectors. The
      79              :  *    reordering rotates the columns of $Z$, keeping them orthonormal and $T$
      80              :  *    upper triangular.
      81              :  * 4. Compute the left eigenvectors as the inverse of the right eigenvectors
      82              :  *    (`dgesv`), so that $\sum_m L^{im} R_{jm} = \delta^{i}_{j}$.
      83              :  *
      84              :  * Finally, we undo the balancing: $R_{im} = D_{mm} R'_{im}$ and
      85              :  * $L^{im} = L'^{im} / D_{mm}$.
      86              :  *
      87              :  * Eigenvalues are treated as equal (degenerate) when they differ by at most
      88              :  * $10^{-10} \|A'\|_F$, where $\|A'\|_F$ is the Frobenius norm of the
      89              :  * balanced matrix.
      90              :  */
      91              : template <size_t MatrixSize>
      92            1 : void find_real_eigensystem(
      93              :     gsl::not_null<tnsr::i<DataVector, MatrixSize>*> eigenvalues,
      94              :     gsl::not_null<tnsr::ij<DataVector, MatrixSize>*> right_eigenvectors,
      95              :     gsl::not_null<tnsr::IJ<DataVector, MatrixSize>*> left_eigenvectors,
      96              :     const tnsr::iJ<DataVector, MatrixSize>& matrix);