Eigenvalues_Module

Classes
struct	Eigen::internal::HessenbergDecompositionMatrixHReturnType< MatrixType >
	Expression type for return value of HessenbergDecomposition::matrixH() More...
struct	Eigen::internal::TridiagonalizationMatrixTReturnType< MatrixType >
	Expression type for return value of Tridiagonalization::matrixT() More...
class	Eigen::ComplexEigenSolver< MatrixType_ >
	Computes eigenvalues and eigenvectors of general complex matrices. More...
class	Eigen::ComplexSchur< MatrixType_ >
	Performs a complex Schur decomposition of a real or complex square matrix. More...
class	Eigen::EigenSolver< MatrixType_ >
	Computes eigenvalues and eigenvectors of general matrices. More...
class	Eigen::GeneralizedEigenSolver< MatrixType_ >
	Computes the generalized eigenvalues and eigenvectors of a pair of general matrices. More...
class	Eigen::GeneralizedSelfAdjointEigenSolver< MatrixType_ >
	Computes eigenvalues and eigenvectors of the generalized selfadjoint eigen problem. More...
class	Eigen::HessenbergDecomposition< MatrixType_ >
	Reduces a square matrix to Hessenberg form by an orthogonal similarity transformation. More...
class	Eigen::RealQZ< MatrixType_ >
	Performs a real QZ decomposition of a pair of square matrices. More...
class	Eigen::RealSchur< MatrixType_ >
	Performs a real Schur decomposition of a square matrix. More...
class	Eigen::SelfAdjointEigenSolver< MatrixType_ >
	Computes eigenvalues and eigenvectors of selfadjoint matrices. More...
class	Eigen::Tridiagonalization< MatrixType_ >
	Tridiagonal decomposition of a selfadjoint matrix. More...

Functions
template<int StorageOrder, typename RealScalar , typename Scalar , typename Index >
static EIGEN_DEVICE_FUNC void	Eigen::internal::tridiagonal_qr_step (RealScalar *diag, RealScalar *subdiag, Index start, Index end, Scalar *matrixQ, Index n)

Detailed Description

Function Documentation

tridiagonal_qr_step()

template<int StorageOrder, typename RealScalar , typename Scalar , typename Index >

static EIGEN_DEVICE_FUNC void Eigen::internal::tridiagonal_qr_step ( RealScalar *  diag, RealScalar *  subdiag, Index  start, Index  end, Scalar *  matrixQ, Index  n  )

static

\eigenvalues_module

Performs a QR step on a tridiagonal symmetric matrix represented as a pair of two vectors diag and subdiag.

Parameters

diag	the diagonal part of the input selfadjoint tridiagonal matrix
subdiag	the sub-diagonal part of the input selfadjoint tridiagonal matrix
start	starting index of the submatrix to work on
end	last+1 index of the submatrix to work on
matrixQ	pointer to the column-major matrix holding the eigenvectors, can be 0
n	size of the input matrix

For compilation efficiency reasons, this procedure does not use eigen expression for its arguments.

Implemented from Golub's "Matrix Computations", algorithm 8.3.2: "implicit symmetric QR step with Wilkinson shift"

                                                                            {
  // Wilkinson Shift.
  RealScalar td = (diag[end - 1] - diag[end]) * RealScalar(0.5);
  RealScalar e = subdiag[end - 1];
  // Note that thanks to scaling, e^2 or td^2 cannot overflow, however they can still
  // underflow thus leading to inf/NaN values when using the following commented code:
  //   RealScalar e2 = numext::abs2(subdiag[end-1]);
  //   RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2));
  // This explain the following, somewhat more complicated, version:
  RealScalar mu = diag[end];
  if (numext::is_exactly_zero(td)) {
    mu -= numext::abs(e);
  } else if (!numext::is_exactly_zero(e)) {
    const RealScalar e2 = numext::abs2(e);
    const RealScalar h = numext::hypot(td, e);
    if (numext::is_exactly_zero(e2)) {
      mu -= e / ((td + (td > RealScalar(0) ? h : -h)) / e);
    } else {
      mu -= e2 / (td + (td > RealScalar(0) ? h : -h));
    }
  }
 
  RealScalar x = diag[start] - mu;
  RealScalar z = subdiag[start];
  // If z ever becomes zero, the Givens rotation will be the identity and
  // z will stay zero for all future iterations.
  for (Index k = start; k < end && !numext::is_exactly_zero(z); ++k) {
    JacobiRotation<RealScalar> rot;
    rot.makeGivens(x, z);
 
    // do T = G' T G
    RealScalar sdk = rot.s() * diag[k] + rot.c() * subdiag[k];
    RealScalar dkp1 = rot.s() * subdiag[k] + rot.c() * diag[k + 1];
 
    diag[k] =
        rot.c() * (rot.c() * diag[k] - rot.s() * subdiag[k]) - rot.s() * (rot.c() * subdiag[k] - rot.s() * diag[k + 1]);
    diag[k + 1] = rot.s() * sdk + rot.c() * dkp1;
    subdiag[k] = rot.c() * sdk - rot.s() * dkp1;
 
    if (k > start) subdiag[k - 1] = rot.c() * subdiag[k - 1] - rot.s() * z;
 
    // "Chasing the bulge" to return to triangular form.
    x = subdiag[k];
    if (k < end - 1) {
      z = -rot.s() * subdiag[k + 1];
      subdiag[k + 1] = rot.c() * subdiag[k + 1];
    }
 
    // apply the givens rotation to the unit matrix Q = Q * G
    if (matrixQ) {
      // FIXME if StorageOrder == RowMajor this operation is not very efficient
      Map<Matrix<Scalar, Dynamic, Dynamic, StorageOrder> > q(matrixQ, n, n);
      q.applyOnTheRight(k, k + 1, rot);
    }
  }
}

References Eigen::numext::abs(), Eigen::numext::abs2(), diag, e(), Eigen::placeholders::end, Eigen::numext::is_exactly_zero(), k, Global_Parameters::mu, n, Eigen::numext::q, rot(), oomph::CumulativeTimings::start(), and plotDoE::x.