/**
- * Return the eigenvalues of a symmetric tensor of rank 2.
+ * Return the eigenvalues of a symmetric 1x1 tensor of rank 2.
+ *
+ * The (single) entry of the tensor is, of course, equal to the (single)
+ * eigenvalue.
*
* @relates SymmetricTensor
* @author Jean-Paul Pelteret, 2017
/**
- * Return the eigenvalues of a symmetric tensor of rank 2.
+ * Return the eigenvalues of a symmetric 2x2 tensor of rank 2.
* The array of eigenvalues is sorted in descending order.
*
* For 2x2 tensors, the eigenvalues of tensor $T$ are the roots of
* <a href="https://en.wikipedia.org/wiki/Eigenvalue_algorithm#2.C3.972_matrices">the characteristic polynomial</a>
- * $0 = \lambda^{2} - \lambda*tr(T) + det(T)$
+ * $0 = \lambda^{2} - \lambda\textrm{tr}(T) + \textrm{det}(T)$
* as given by
- * $\lambda = \frac{tr(T) \pm \sqrt{tr^{2}(T) - 4*det(T)}}{2}$.
+ * $\lambda = \frac{\textrm{tr}(T) \pm \sqrt{[\textrm{tr}(T)]^{2} - 4\textrm{det}(T)}}{2}$.
*
* @warning The algorithm employed here determines the eigenvalues by
* computing the roots of the characteristic polynomial. In the case that there
* exists a common root (the eigenvalues are equal), the computation is
* <a href="https://scicomp.stackexchange.com/q/23686">subject to round-off errors</a>
* of order $\sqrt{\epsilon}$.
+ * As an alternative, the eigenvectors() function provides a more robust, but costly,
+ * method to compute the eigenvalues of a symmetric tensor.
*
* @relates SymmetricTensor
* @author Jean-Paul Pelteret, 2017
/**
- * Return the eigenvalues of a symmetric tensor of rank 2.
+ * Return the eigenvalues of a symmetric 3x3 tensor of rank 2.
* The array of eigenvalues is sorted in descending order.
*
* For 3x3 tensors, the eigenvalues of tensor $T$ are the roots of
* <a href="https://en.wikipedia.org/wiki/Eigenvalue_algorithm#3.C3.973_matrices">the characteristic polynomial</a>
- * $0 = \lambda^{3} - \lambda^{2}*tr(T) - \frac{1}{2} \lambda (tr(T^{2}) - tr^{2}(T)) - det(T)$.
+ * $0 = \lambda^{3} - \lambda^{2}\textrm{tr}(T) - \frac{1}{2} \lambda (\textrm{tr}(T^{2}) - [\textrm{tr}(T)]^{2}) - \textrm{det}(T)$.
*
* @warning The algorithm employed here determines the eigenvalues by
* computing the roots of the characteristic polynomial. In the case that there
* exists a common root (the eigenvalues are equal), the computation is
* <a href="https://scicomp.stackexchange.com/q/23686">subject to round-off errors</a>
* of order $\sqrt{\epsilon}$.
+ * As an alternative, the eigenvectors() function provides a more robust, but costly,
+ * method to compute the eigenvalues of a symmetric tensor.
*
* @relates SymmetricTensor
* @author Jean-Paul Pelteret, 2017
+namespace internal
+{
+ namespace
+ {
+ /**
+ * Tridiagonalize a rank-2 symmetric using the Householder method.
+ * The specialized algorithm implemented here is given in
+ * Kopp, J.
+ * Efficient numerical diagonalization of hermitian 3x3 matrices
+ * International Journal of Modern Physics C, 2008, 19, 523-548
+ * doi: 10.1142/S0129183108012303
+ * arXiv.org preprint: physics/0610206
+ * and is based off of the generic algorithm presented in section 11.3.2 of
+ * Press, W. H.
+ * Numerical recipes 3rd edition: The art of scientific computing
+ * Cambridge university press, 2007
+ *
+ * @param[in] A This tensor to be tridiagonalized
+ * @param[out] Q The orthogonal matrix effecting the transformation
+ * @param[out] d The diagonal elements of the tridiagonal matrix
+ * @param[out] e The off-diagonal elements of the tridiagonal matrix
+ */
+ template <int dim, typename Number>
+ void
+ tridiagonalize (const dealii::SymmetricTensor<2,dim,Number> &A,
+ dealii::Tensor<2,dim,Number> &Q,
+ std::array<Number,dim> &d,
+ std::array<Number,dim-1> &e)
+ {
+ // Create some intermediate storage
+ Number h,g,omega_inv,K,f;
+
+ // Initialize the transformation matrix as the
+ // identity tensor
+ Q = dealii::unit_symmetric_tensor<dim,Number>();
+
+ // Make the first row and column to be of the
+ // desired form
+ h = 0.0;
+ for (int i=1; i < dim; i++)
+ h += A[0][i]*A[0][i];
+
+ g = 0.0;
+ if (A[0][1] > 0.0)
+ g = -std::sqrt(h);
+ else
+ g = std::sqrt(h);
+ e[0] = g;
+
+ std::array<Number,dim> u;
+ for (int i=1; i < dim; i++)
+ {
+ u[i] = A[0][i];
+ if (i == 1)
+ u[i] -= g;
+ }
+
+ std::array<Number,dim> q;
+ const Number omega = h - g * A[0][1];
+ if (omega > 0.0)
+ {
+ omega_inv = 1.0 / omega;
+ K = 0.0;
+ for (int i=1; i < dim; i++)
+ {
+ f = 0.0;
+ for (int j=1; j < dim; j++)
+ f += A[i][j] * u[j];
+ q[i] = omega_inv * f;
+ K += u[i] * f;
+ }
+ K *= 0.5*omega_inv*omega_inv;
+
+ for (int i=1; i < dim; i++)
+ q[i] = q[i] - K * u[i];
+
+ d[0] = A[0][0];
+ for (int i=1; i < dim; i++)
+ d[i] = A[i][i] - 2.0*q[i]*u[i];
+
+ // Store inverse Householder transformation
+ // in Q
+ for (int j=1; j < dim; j++)
+ {
+ f = omega_inv * u[j];
+ for (int i=1; i < dim; i++)
+ Q[i][j] = Q[i][j] - f*u[i];
+ }
+
+ // For dim = 3: Calculate updated A[1][2] and
+ // store it in e[1]
+ for (int i=1; i < dim-1; i++)
+ e[i] = A[i][i+1] - q[i]*u[i+1] - u[i]*q[i+1];
+ }
+ else
+ {
+ for (int i=0; i < dim; i++)
+ d[i] = A[i][i];
+
+ // For dim = 3:
+ for (int i=1; i < dim-1; i++)
+ e[i] = A[i][i+1];
+ }
+ }
+
+
+ /**
+ * Compute the eigenvalues and eigenvectors of a real-valued rank-2
+ * symmetric tensor using the QL algorithm with implicit shifts.
+ * The specialized algorithm implemented here is given in
+ * Kopp, J.
+ * Efficient numerical diagonalization of hermitian 3x3 matrices
+ * International Journal of Modern Physics C, 2008, 19, 523-548
+ * doi: 10.1142/S0129183108012303
+ * arXiv.org preprint: physics/0610206
+ * and is based off of the generic algorithm presented in section 11.4.3 of
+ * Press, W. H.
+ * Numerical recipes 3rd edition: The art of scientific computing
+ * Cambridge university press, 2007.
+ *
+ * @param[in] A The tensor of which the eigenvectors and eigenvalues are
+ * to be computed.
+ *
+ * @return An array containing the eigenvectors and the associated eigenvalues
+ */
+ template <int dim, typename Number>
+ std::array<std::pair<Number, Tensor<1,dim,Number> >,dim>
+ ql_implicit_shifts (const dealii::SymmetricTensor<2,dim,Number> &A)
+ {
+ static_assert(numbers::NumberTraits<Number>::is_complex == false,
+ "This implementation of the QL implicit shift algorithm does "
+ "not support complex numbers");
+
+ // Transform A to real tridiagonal form by the Householder method:
+ // The orthogonal matrix effecting the transformation
+ // this will ultimately store the eigenvectors
+ dealii::Tensor<2,dim,Number> Q;
+ // The diagonal elements of the tridiagonal matrix;
+ // this will ultimately store the eigenvalues
+ std::array<Number,dim> w;
+ // The off-diagonal elements of the tridiagonal
+ std::array<Number,dim-1> ee;
+ tridiagonalize(A, Q, w, ee);
+
+ // Number of iterations
+ const unsigned int max_n_it = 30;
+
+ // Transfer the off-diagonal entries to an auxiliary array
+ // The third element is used only as temporary workspace
+ std::array<Number,dim> e;
+ for (unsigned int i=0; i<dim-1; ++i)
+ e[i] = ee[i];
+
+ // Create some intermediate storage
+ Number g, r, p, f, b, s, c, t;
+
+ // Loop over all off-diagonal elements
+ for (int l=0; l < dim-1; l++)
+ {
+ for (unsigned int it=0; it <= max_n_it; ++it)
+ {
+ // Check for convergence and exit iteration loop
+ // if the off-diagonal element e[l] is zero
+ int m = l;
+ for (; m <= dim-2; m++)
+ {
+ g = std::abs(w[m]) + std::abs(w[m+1]);
+ if (std::abs(e[m]) + g == g)
+ break;
+ }
+ if (m == l)
+ break;
+
+ // Throw if no convergence is achieved within a
+ // stipulated number of iterations
+ if (it == max_n_it)
+ {
+ AssertThrow(false, ExcMessage("No convergence in iterative QL eigenvector algorithm."))
+ return std::array<std::pair<Number, Tensor<1,dim,Number> >,dim> ();
+ }
+
+ // Calculate the shift..
+ g = (w[l+1] - w[l]) / (e[l] + e[l]);
+ r = std::sqrt(g*g + 1.0);
+ // .. and then compute g = d_m - k_s for the
+ // plane rotation (Press2007a eq 11.4.22)
+ if (g > 0.0)
+ g = w[m] - w[l] + e[l]/(g + r);
+ else
+ g = w[m] - w[l] + e[l]/(g - r);
+
+ // Perform plane rotation, as is done in the
+ // standard QL algorithm, followed by Givens
+ // rotations to recover the tridiagonal form
+ s = c = 1.0;
+ p = 0.0;
+ for (int i=m-1; i >= l; i--)
+ {
+ f = s * e[i];
+ b = c * e[i];
+
+ // Branch to recover from underflow
+ if (std::abs(f) > std::abs(g))
+ {
+ c = g / f;
+ r = std::sqrt(c*c + 1.0);
+ e[i+1] = f * r;
+ c *= (s = 1.0/r);
+ }
+ else
+ {
+ s = f / g;
+ r = std::sqrt(s*s + 1.0);
+ e[i+1] = g * r;
+ s *= (c = 1.0/r);
+ }
+
+ g = w[i+1] - p;
+ r = (w[i] - g)*s + 2.0*c*b;
+ p = s * r;
+ w[i+1] = g + p;
+ g = c*r - b;
+
+ // Form the eigenvectors
+ for (int k=0; k < dim; k++)
+ {
+ t = Q[k][i+1];
+ Q[k][i+1] = s*Q[k][i] + c*t;
+ Q[k][i] = c*Q[k][i] - s*t;
+ }
+ }
+ w[l] -= p;
+ e[l] = g;
+ e[m] = 0.0;
+ }
+ }
+
+ // Structure the data to be outputted
+ std::array<std::pair<Number, Tensor<1,dim,Number> >,dim> eig_vals_vecs;
+ for (unsigned int e=0; e<dim; ++e)
+ {
+ eig_vals_vecs[e].first = w[e];
+
+ // The column "e" of Q contains the non-normalized
+ // eigenvector associated with the eigenvalue "e"
+ for (unsigned int a=0; a<dim; ++a)
+ {
+ eig_vals_vecs[e].second[a] = Q[a][e];
+ }
+
+ // Normalize
+ Assert(eig_vals_vecs[e].second.norm() != 0.0, ExcDivideByZero());
+ eig_vals_vecs[e].second /= eig_vals_vecs[e].second.norm();
+ }
+ return eig_vals_vecs;
+ }
+
+
+ /**
+ * Compute the eigenvalues and eigenvectors of a real-valued rank-2
+ * symmetric tensor using the Jacobi algorithm.
+ * The specialized algorithm implemented here is given in
+ * Kopp, J.
+ * Efficient numerical diagonalization of hermitian 3x3 matrices
+ * International Journal of Modern Physics C, 2008, 19, 523-548
+ * doi: 10.1142/S0129183108012303
+ * arXiv.org preprint: physics/0610206
+ * and is based off of the generic algorithm presented in section 11.4.3 of
+ * Press, W. H.
+ * Numerical recipes 3rd edition: The art of scientific computing
+ * Cambridge university press, 2007
+ *
+ * @param[in] A The tensor of which the eigenvectors and eigenvalues are
+ * to be computed.
+ *
+ * @return An array containing the eigenvectors and the associated eigenvalues
+ */
+ template <int dim, typename Number>
+ std::array<std::pair<Number, Tensor<1,dim,Number> >,dim>
+ jacobi (dealii::SymmetricTensor<2,dim,Number> A)
+ {
+ static_assert(numbers::NumberTraits<Number>::is_complex == false,
+ "This implementation of the Jacobi algorithm does "
+ "not support complex numbers");
+
+ // Sums of diagonal resp. off-diagonal elements
+ Number sd, so;
+ // sin(phi), cos(phi), tan(phi) and temporary storage
+ Number s, c, t;
+ // More temporary storage
+ Number g, h, z, theta;
+ // Threshold value
+ Number thresh;
+
+ // Initialize the transformation matrix as the
+ // identity tensor
+ dealii::Tensor<2,dim,Number> Q (dealii::unit_symmetric_tensor<dim,Number>());
+
+ // The diagonal elements of the tridiagonal matrix;
+ // this will ultimately store the eigenvalues
+ std::array<Number,dim> w;
+ for (int i=0; i < dim; i++)
+ w[i] = A[i][i];
+
+ // Calculate (tr(A))^{2}
+ sd = trace(A);
+ sd *= sd;
+
+ // Number of iterations
+ const unsigned int max_n_it = 50;
+ for (unsigned int it=0; it <= max_n_it; it++)
+ {
+ // Test for convergence
+ so = 0.0;
+ for (int p=0; p < dim; p++)
+ for (int q=p+1; q < dim; q++)
+ so += std::abs(A[p][q]);
+ if (so == 0.0)
+ break;
+
+ // Throw if no convergence is achieved within a
+ // stipulated number of iterations
+ if (it == max_n_it)
+ {
+ AssertThrow(false, ExcMessage("No convergence in iterative Jacobi eigenvector algorithm."))
+ return std::array<std::pair<Number, Tensor<1,dim,Number> >,dim> ();
+ }
+
+ // Compute threshold value which dictates whether or
+ // not a Jacobi rotation is performed
+ const unsigned int n_it_skip = 4;
+ if (it < n_it_skip)
+ thresh = 0.2 * so / (dim*dim);
+ else
+ thresh = 0.0;
+
+ // Perform sweep
+ for (int p=0; p < dim; p++)
+ for (int q=p+1; q < dim; q++)
+ {
+ g = 100.0 * std::abs(A[p][q]);
+
+ // After a given number of iterations the
+ // rotation is skipped if the off-diagonal
+ // element is small
+ if (it > n_it_skip &&
+ std::abs(w[p]) + g == std::abs(w[p]) &&
+ std::abs(w[q]) + g == std::abs(w[q]))
+ {
+ A[p][q] = 0.0;
+ }
+ else if (std::abs(A[p][q]) > thresh)
+ {
+ // Calculate Jacobi transformation
+ h = w[q] - w[p];
+
+ // Compute surrogate for angle theta resulting from
+ // angle transformation and subsequent smallest solution
+ // of quadratic equation
+ if (std::abs(h) + g == std::abs(h))
+ {
+ // Prevent overflow for large theta^2. This computation
+ // is the algebraic equivalent of t = 1/(2*theta).
+ t = A[p][q] / h;
+ }
+ else
+ {
+ theta = 0.5 * h / A[p][q];
+ if (theta < 0.0)
+ t = -1.0 / (std::sqrt(1.0 + theta*theta) - theta);
+ else
+ t = 1.0 / (std::sqrt(1.0 + theta*theta) + theta);
+ }
+
+ // Compute trigonometric functions for rotation
+ // in such a way as to prevent overflow for
+ // large theta.
+ c = 1.0/std::sqrt(1.0 + t*t);
+ s = t * c;
+ z = t * A[p][q];
+
+ // Apply Jacobi transformation...
+ A[p][q] = 0.0;
+ w[p] -= z;
+ w[q] += z;
+ // ... by executing the various rotations in sequence
+ for (int r=0; r < p; r++)
+ {
+ t = A[r][p];
+ A[r][p] = c*t - s*A[r][q];
+ A[r][q] = s*t + c*A[r][q];
+ }
+ for (int r=p+1; r < q; r++)
+ {
+ t = A[p][r];
+ A[p][r] = c*t - s*A[r][q];
+ A[r][q] = s*t + c*A[r][q];
+ }
+ for (int r=q+1; r < dim; r++)
+ {
+ t = A[p][r];
+ A[p][r] = c*t - s*A[q][r];
+ A[q][r] = s*t + c*A[q][r];
+ }
+
+ // Update the eigenvectors
+ for (int r=0; r < dim; r++)
+ {
+ t = Q[r][p];
+ Q[r][p] = c*t - s*Q[r][q];
+ Q[r][q] = s*t + c*Q[r][q];
+ }
+ }
+ }
+ }
+
+ // Structure the data to be outputted
+ std::array<std::pair<Number, Tensor<1,dim,Number> >,dim> eig_vals_vecs;
+ for (unsigned int e=0; e<dim; ++e)
+ {
+ eig_vals_vecs[e].first = w[e];
+
+ // The column "e" of Q contains the non-normalized
+ // eigenvector associated with the eigenvalue "e"
+ for (unsigned int a=0; a<dim; ++a)
+ {
+ eig_vals_vecs[e].second[a] = Q[a][e];
+ }
+
+ // Normalize
+ Assert(eig_vals_vecs[e].second.norm() != 0.0, ExcDivideByZero());
+ eig_vals_vecs[e].second /= eig_vals_vecs[e].second.norm();
+ }
+ return eig_vals_vecs;
+ }
+
+
+
+ /**
+ * Compute the eigenvalues and eigenvectors of a real-valued rank-2
+ * symmetric tensor using the characteristic equation to compute eigenvalues
+ * and an analytical approach based on the cross-product for the eigenvectors.
+ * If the computations are deemed too inaccurate then the method falls
+ * back to ql_implicit_shifts.
+ *
+ * @param[in] A The tensor of which the eigenvectors and eigenvalues are
+ * to be computed.
+ *
+ * @return An array containing the eigenvectors and the associated eigenvalues
+ */
+ template <typename Number>
+ std::array<std::pair<Number, Tensor<1,2,Number> >,2>
+ hybrid (const dealii::SymmetricTensor<2,2,Number> &A)
+ {
+ static_assert(numbers::NumberTraits<Number>::is_complex == false,
+ "This implementation of the 2d Hybrid algorithm does "
+ "not support complex numbers");
+
+ const unsigned int dim = 2;
+
+ // Calculate eigenvalues
+ const std::array<Number,dim> w = eigenvalues(A);
+
+ std::array<std::pair<Number, Tensor<1,dim,Number> >,dim> eig_vals_vecs;
+
+ Number t, u; // Intermediate storage
+ t = std::abs(w[0]);
+ for (unsigned int i=1; i<dim; ++i)
+ {
+ u = std::abs(w[i]);
+ if (u > t)
+ t = u;
+ }
+
+ if (t < 1.0)
+ u = t;
+ else
+ u = t*t;
+
+ // Estimated maximum roundoff error
+ const double error = 256.0 * std::numeric_limits<double>::epsilon() * u*u;
+
+ // Store eigenvalues
+ eig_vals_vecs[0].first = w[0];
+ eig_vals_vecs[1].first = w[1];
+
+ // Compute eigenvectors
+ // http://www.math.harvard.edu/archive/21b_fall_04/exhibits/2dmatrices/
+ // https://math.stackexchange.com/a/1548616
+ if (A[1][0] != 0.0)
+ {
+ // First eigenvector
+ eig_vals_vecs[0].second[0] = w[0] - A[1][1];
+ eig_vals_vecs[0].second[1] = A[1][0];
+
+ // Second eigenvector
+ eig_vals_vecs[1].second[0] = w[1] - A[1][1];
+ eig_vals_vecs[1].second[1] = A[1][0];
+ }
+ else
+ {
+ // First eigenvector
+ eig_vals_vecs[0].second[0] = w[0];
+ eig_vals_vecs[0].second[1] = 0.0;
+
+ // Second eigenvector
+ eig_vals_vecs[1].second[0] = 0.0;
+ eig_vals_vecs[1].second[1] = w[1];
+ }
+ // Normalize
+ eig_vals_vecs[0].second /= eig_vals_vecs[0].second.norm();
+ eig_vals_vecs[1].second /= eig_vals_vecs[1].second.norm();
+
+ // If vectors are nearly linearly dependent, or if there might have
+ // been large cancelations in the calculation of A[i][i] - w[0], fall
+ // back to QL algorithm
+ if (eig_vals_vecs[0].second * eig_vals_vecs[1].second > error)
+ {
+ return ql_implicit_shifts(A);
+ }
+
+ return eig_vals_vecs;
+ }
+
+
+
+ /**
+ * Compute the eigenvalues and eigenvectors of a real-valued rank-2
+ * symmetric tensor using the characteristic equation to compute eigenvalues
+ * and an analytical approach based on the cross-product for the eigenvectors.
+ * If the computations are deemed too inaccurate then the method falls
+ * back to ql_implicit_shifts.
+ * The specialized algorithm implemented here is given in
+ * Kopp, J.
+ * Efficient numerical diagonalization of hermitian 3x3 matrices
+ * International Journal of Modern Physics C, 2008, 19, 523-548
+ * doi: 10.1142/S0129183108012303
+ * arXiv.org preprint: physics/0610206
+ *
+ * @param[in] A The tensor of which the eigenvectors and eigenvalues are
+ * to be computed.
+ *
+ * @return An array containing the eigenvectors and the associated eigenvalues
+ */
+ template <typename Number>
+ std::array<std::pair<Number, Tensor<1,3,Number> >,3>
+ hybrid (const dealii::SymmetricTensor<2,3,Number> &A)
+ {
+ static_assert(numbers::NumberTraits<Number>::is_complex == false,
+ "This implementation of the 3d Hybrid algorithm does "
+ "not support complex numbers");
+
+ const unsigned int dim = 3;
+ Number norm; // Squared norm or inverse norm of current eigenvector
+ Number t, u; // Intermediate storage
+
+ // Calculate eigenvalues
+ const std::array<Number,dim> w = eigenvalues(A);
+
+ t = std::abs(w[0]);
+ for (unsigned int i=1; i<dim; ++i)
+ {
+ u = std::abs(w[i]);
+ if (u > t)
+ t = u;
+ }
+
+ if (t < 1.0)
+ u = t;
+ else
+ u = t*t;
+
+ // Estimated maximum roundoff error
+ const double error = 256.0 * std::numeric_limits<double>::epsilon() * u*u;
+
+ // Initialize the transformation matrix as the
+ // identity tensor
+ dealii::Tensor<2,dim,Number> Q;
+ Q[0][1] = A[0][1]*A[1][2] - A[0][2]*A[1][1];
+ Q[1][1] = A[0][2]*A[0][1] - A[1][2]*A[0][0];
+ Q[2][1] = A[0][1]*A[0][1];
+
+ // Calculate first eigenvector by the formula
+ // v[0] = (A - w[0]).e1 x (A - w[0]).e2
+ Q[0][0] = Q[0][1] + A[0][2]*w[0];
+ Q[1][0] = Q[1][1] + A[1][2]*w[0];
+ Q[2][0] = (A[0][0] - w[0]) * (A[1][1] - w[0]) - Q[2][1];
+ norm = Q[0][0]*Q[0][0] + Q[1][0]*Q[1][0] + Q[2][0]*Q[2][0];
+
+ // If vectors are nearly linearly dependent, or if there might have
+ // been large cancellations in the calculation of A[i][i] - w[0], fall
+ // back to QL algorithm
+ // Note that this simultaneously ensures that multiple eigenvalues do
+ // not cause problems: If w[0] = w[1], then A - w[0] * I has rank 1,
+ // i.e. all columns of A - w[0] * I are linearly dependent.
+ if (norm <= error)
+ {
+ return ql_implicit_shifts(A);
+ }
+ else // This is the standard branch
+ {
+ norm = std::sqrt(1.0 / norm);
+ for (unsigned j=0; j < dim; j++)
+ Q[j][0] = Q[j][0] * norm;
+ }
+
+ // Calculate second eigenvector by the formula
+ // v[1] = (A - w[1]).e1 x (A - w[1]).e2
+ Q[0][1] = Q[0][1] + A[0][2]*w[1];
+ Q[1][1] = Q[1][1] + A[1][2]*w[1];
+ Q[2][1] = (A[0][0] - w[1]) * (A[1][1] - w[1]) - Q[2][1];
+ norm = Q[0][1]*Q[0][1] + Q[1][1]*Q[1][1] + Q[2][1]*Q[2][1];
+ if (norm <= error)
+ {
+ return ql_implicit_shifts(A);
+ }
+ else
+ {
+ norm = std::sqrt(1.0 / norm);
+ for (unsigned int j=0; j < dim; j++)
+ Q[j][1] = Q[j][1] * norm;
+ }
+
+ // Calculate third eigenvector according to
+ // v[2] = v[0] x v[1]
+ Q[0][2] = Q[1][0]*Q[2][1] - Q[2][0]*Q[1][1];
+ Q[1][2] = Q[2][0]*Q[0][1] - Q[0][0]*Q[2][1];
+ Q[2][2] = Q[0][0]*Q[1][1] - Q[1][0]*Q[0][1];
+
+ // Structure the data to be outputted
+ std::array<std::pair<Number, Tensor<1,dim,Number> >,dim> eig_vals_vecs;
+ for (unsigned int e=0; e<dim; ++e)
+ {
+ eig_vals_vecs[e].first = w[e];
+
+ // The column "e" of Q contains the non-normalized
+ // eigenvector associated with the eigenvalue "e"
+ for (unsigned int a=0; a<dim; ++a)
+ {
+ eig_vals_vecs[e].second[a] = Q[a][e];
+ }
+
+ // Normalize
+ Assert(eig_vals_vecs[e].second.norm() != 0.0, ExcDivideByZero());
+ eig_vals_vecs[e].second /= eig_vals_vecs[e].second.norm();
+ }
+ return eig_vals_vecs;
+ }
+
+
+ template<int dim, typename Number>
+ struct SortEigenValuesVectors
+ {
+ typedef std::pair<Number, Tensor<1,dim,Number> > EigValsVecs;
+ bool operator() (const EigValsVecs &lhs,
+ const EigValsVecs &rhs)
+ {
+ return lhs.first > rhs.first;
+ }
+ };
+
+ }
+} // namespace internal
+
+
+
+/**
+ * An enumeration for the algorithm to be employed when performing
+ * the computation of normalized eigenvectors and their corresponding
+ * eigenvalues.
+ *
+ * The specialized algorithms utilized in computing the eigenvectors are
+ * presented in
+ * @code{.bib}
+ * @Article{Kopp2008,
+ * title = {Efficient numerical diagonalization of hermitian 3x3 matrices},
+ * author = {Kopp, J.},
+ * journal = {International Journal of Modern Physics C},
+ * year = {2008},
+ * volume = {19},
+ * number = {3},
+ * pages = {523--548},
+ * doi = {10.1142/S0129183108012303},
+ * eprinttype = {arXiv},
+ * eprint = {physics/0610206v3},
+ * eprintclass = {physics.comp-ph},
+ * url = {https://www.mpi-hd.mpg.de/personalhomes/globes/3x3/index.html}
+ * }
+ * @endcode
+ */
+enum EigenvectorMethod
+{
+ /**
+ * A hybrid approach that preferentially uses the characteristic equation to
+ * compute eigenvalues and an analytical approach based on the cross-product
+ * for the eigenvectors. If the computations are deemed too inaccurate then
+ * the method falls back to ql_implicit_shifts.
+ *
+ * This method potentially offers the quickest computation if the pathological
+ * case is not encountered.
+ */
+ hybrid,
+ /**
+ * The iterative QL algorithm with implicit shifts applied after
+ * tridiagonalisation of the tensor using the householder method.
+ *
+ * This method offers a compromise between speed of computation and its
+ * robustness. This method is particularly useful when the elements
+ * of $T$ have greatly varying magnitudes, which would typically lead to a
+ * loss of accuracy when computing the smaller eigenvalues.
+ */
+ ql_implicit_shifts,
+ /**
+ * The iterative Jacobi algorithm.
+ *
+ * This method offers is the most robust of the available options, with
+ * reliable results obtained for even the most pathological cases. It is,
+ * however, the slowest algorithm of all of those implemented.
+ */
+ jacobi
+};
+
+
+
+/**
+ * Return the eigenvalues and eigenvectors of a symmetric tensor of rank 2.
+ *
+ * @relates SymmetricTensor
+ * @author Jean-Paul Pelteret, 2017
+ */
+template <typename Number>
+std::array<std::pair<Number, Tensor<1,1,Number> >,1>
+eigenvectors (const SymmetricTensor<2,1,Number> &T,
+ const enum EigenvectorMethod /*method*/)
+{
+ return { {std::make_pair(T[0][0], Tensor<1,1,Number>({1.0}))} };
+}
+
+
+
+/**
+ * Return the eigenvalues and eigenvectors of a real-valued rank-2 symmetric
+ * tensor $T$.
+ *
+ * The specialized algorithms utilized in computing the eigenvectors are
+ * presented in
+ * @code{.bib}
+ * @Article{Kopp2008,
+ * title = {Efficient numerical diagonalization of hermitian 3x3 matrices},
+ * author = {Kopp, J.},
+ * journal = {International Journal of Modern Physics C},
+ * year = {2008},
+ * volume = {19},
+ * number = {3},
+ * pages = {523--548},
+ * doi = {10.1142/S0129183108012303},
+ * eprinttype = {arXiv},
+ * eprint = {physics/0610206v3},
+ * eprintclass = {physics.comp-ph},
+ * url = {https://www.mpi-hd.mpg.de/personalhomes/globes/3x3/index.html}
+ * }
+ * @endcode
+ *
+ * @relates SymmetricTensor
+ * @author Joachim Kopp, Jean-Paul Pelteret, 2017
+ */
+template <typename Number>
+std::array<std::pair<Number, Tensor<1,2,Number> >,2>
+eigenvectors (const SymmetricTensor<2,2,Number> &T,
+ const enum EigenvectorMethod method = ql_implicit_shifts)
+{
+ std::array<std::pair<Number, Tensor<1,2,Number> >,2> eig_vals_vecs;
+
+ if (method == hybrid)
+ eig_vals_vecs = internal::hybrid(T);
+ else if (method == ql_implicit_shifts)
+ eig_vals_vecs = internal::ql_implicit_shifts(T);
+ else if (method == jacobi)
+ eig_vals_vecs = internal::jacobi(T);
+ else
+ AssertThrow(false, ExcNotImplemented());
+
+ std::sort(eig_vals_vecs.begin(), eig_vals_vecs.end(),
+ internal::SortEigenValuesVectors<2,Number>());
+ return eig_vals_vecs;
+}
+
+
+
+/**
+ * Return the eigenvalues and eigenvectors of a real-valued rank-2 symmetric
+ * tensor $T$.
+ *
+ * The specialized algorithms utilized in computing the eigenvectors are
+ * presented in
+ * @code{.bib}
+ * @Article{Kopp2008,
+ * title = {Efficient numerical diagonalization of hermitian 3x3 matrices},
+ * author = {Kopp, J.},
+ * journal = {International Journal of Modern Physics C},
+ * year = {2008},
+ * volume = {19},
+ * number = {3},
+ * pages = {523--548},
+ * doi = {10.1142/S0129183108012303},
+ * eprinttype = {arXiv},
+ * eprint = {physics/0610206v3},
+ * eprintclass = {physics.comp-ph},
+ * url = {https://www.mpi-hd.mpg.de/personalhomes/globes/3x3/index.html}
+ * }
+ * @endcode
+ *
+ * @relates SymmetricTensor
+ * @author Joachim Kopp, Jean-Paul Pelteret, 2017
+ */
+template <typename Number>
+std::array<std::pair<Number, Tensor<1,3,Number> >,3>
+eigenvectors (const SymmetricTensor<2,3,Number> &T,
+ const enum EigenvectorMethod method = ql_implicit_shifts)
+{
+ std::array<std::pair<Number, Tensor<1,3,Number> >,3> eig_vals_vecs;
+
+ if (method == hybrid)
+ eig_vals_vecs = internal::hybrid(T);
+ else if (method == ql_implicit_shifts)
+ eig_vals_vecs = internal::ql_implicit_shifts(T);
+ else if (method == jacobi)
+ eig_vals_vecs = internal::jacobi(T);
+ else
+ AssertThrow(false, ExcNotImplemented());
+
+ std::sort(eig_vals_vecs.begin(), eig_vals_vecs.end(),
+ internal::SortEigenValuesVectors<3,Number>());
+ return eig_vals_vecs;
+}
+
+
+
/**
* Return the transpose of the given symmetric tensor. Since we are working
* with symmetric objects, the transpose is of course the same as the original
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Test that SymmetricTensor eigenector calculations are correct
+
+#include "../tests.h"
+
+#include <deal.II/base/symmetric_tensor.h>
+#include <array>
+
+void
+check_value (const int dim, const int index,
+ const double expected, const double actual,
+ const double tol = 1e-12)
+{
+ const double rel_error = std::abs(expected - actual)/std::abs(actual);
+ if (rel_error > tol)
+ {
+ deallog
+ << "Incorrect eigenvalue calculated: "
+ << "Dim " << dim
+ << ", index " << index
+ << ". Expected " << expected
+ << ", actual: " << actual
+ << ", relative error: " << rel_error
+ << std::endl;
+ }
+};
+
+template<int dim>
+bool
+is_unit_vector(const Tensor<1,dim> &v)
+{
+ return std::abs(v.norm() - 1.0) < 1e-9;
+}
+
+template<int dim>
+bool
+check_orientation(Tensor<1,dim> v1,
+ Tensor<1,dim> v2,
+ const double tol = 1e-9)
+{
+ v1 /= v1.norm();
+ v2 /= v2.norm();
+ return std::abs(std::abs(v1*v2) - 1.0) < tol;
+}
+
+template<int dim>
+void
+check_vector (const int index,
+ const Tensor<1,dim> expected, const Tensor<1,dim> actual,
+ const double tol = 1e-12)
+{
+ const bool orientation = check_orientation(expected,actual);
+ const bool unit_vec = is_unit_vector(actual);
+ if (!(orientation & unit_vec))
+ {
+ deallog
+ << "Incorrect eigenvector calculated: "
+ << "Dim " << dim
+ << ", index " << index
+ << ". Expected " << expected
+ << ", actual: " << actual
+ << ", orientation check: " << orientation
+ << ", unit vector check: " << unit_vec
+ << std::endl;
+ }
+};
+
+void
+test_dim_1 (const enum EigenvectorMethod method,
+ const double e1, const double tol = 1e-12)
+{
+ const unsigned int dim = 1;
+ SymmetricTensor<2,dim> T;
+ T[0][0] = e1;
+ const auto eig_vals_vecs = eigenvectors(T, method);
+
+ check_value(dim,0, e1, eig_vals_vecs[0].first, tol);
+ check_vector(0, Tensor<1,dim>({1}), eig_vals_vecs[0].second);
+}
+
+void
+test_dim_2 (const enum EigenvectorMethod method,
+ const double e1, Tensor<1,2> v1,
+ const double e2, const double tol = 1e-12)
+{
+ const unsigned int dim = 2;
+ v1 /= v1.norm();
+ const Tensor<1,dim> v2 = cross_product_2d(v1);
+
+ Assert(is_unit_vector(v1), ExcMessage("Vector is not of unit length."));
+ Assert(is_unit_vector(v2), ExcMessage("Vector is not of unit length."));
+ Assert(e1 >= e2, ExcMessage("Input eigenvalue ordering is not correct."));
+
+ const SymmetricTensor<2,dim> T
+ = e1*symmetrize(outer_product(v1,v1))
+ + e2*symmetrize(outer_product(v2,v2));
+
+ const auto eig_vals_vecs = eigenvectors(T, method);
+
+ check_value(dim,0, e1, eig_vals_vecs[0].first, tol);
+ check_value(dim,1, e2, eig_vals_vecs[1].first, tol);
+ check_vector(0, v1, eig_vals_vecs[0].second);
+ check_vector(1, v2, eig_vals_vecs[1].second);
+}
+
+void
+test_dim_3 (const enum EigenvectorMethod method,
+ const double e1, Tensor<1,3> v1,
+ const double e2, Tensor<1,3> v2,
+ const double e3, const double tol = 1e-12)
+{
+ const unsigned int dim = 3;
+
+ // Note: We do not necessarily expect the seed directors v1,v2 to be orthogonal
+ v1 /= v1.norm();
+ v2 /= v2.norm();
+ Tensor<1,dim> v3 = cross_product_3d(v1,v2);
+ v3 /= v3.norm();
+ v1 = cross_product_3d(v2,v3);
+
+ Assert(is_unit_vector(v1), ExcMessage("Vector is not of unit length."));
+ Assert(is_unit_vector(v2), ExcMessage("Vector is not of unit length."));
+ Assert(is_unit_vector(v3), ExcMessage("Vector is not of unit length."));
+ Assert(check_orientation(v2, cross_product_3d(v3,v1)), ExcMessage("Vectors are not orthogonal."));
+ Assert(e1 >= e2, ExcMessage("Input eigenvalue ordering is not correct."));
+ Assert(e2 >= e3, ExcMessage("Input eigenvalue ordering is not correct."));
+
+ const SymmetricTensor<2,dim> T
+ = e1*symmetrize(outer_product(v1,v1))
+ + e2*symmetrize(outer_product(v2,v2))
+ + e3*symmetrize(outer_product(v3,v3));
+
+ const auto eig_vals_vecs = eigenvectors(T, method);
+
+ check_value(dim,0, e1, eig_vals_vecs[0].first, tol);
+ check_value(dim,1, e2, eig_vals_vecs[1].first, tol);
+ check_value(dim,2, e3, eig_vals_vecs[2].first, tol);
+ check_vector(0, v1, eig_vals_vecs[0].second);
+ check_vector(1, v2, eig_vals_vecs[1].second);
+ check_vector(2, v3, eig_vals_vecs[2].second);
+}
+
+
+void run_tests(const enum EigenvectorMethod method)
+{
+ // Dim = 1
+ {
+ deallog.push("Test 1");
+ {
+ test_dim_1(method, 3.6);
+ }
+ deallog.pop();
+ }
+
+ // Dim = 2
+ {
+ // Diagonal
+ deallog.push("Test 2a");
+ {
+ test_dim_2(method,
+ 3.6, Tensor<1,2>({1,0}),
+ 2.4 );
+ }
+ deallog.pop();
+
+ // Diagonal (large difference)
+ deallog.push("Test 2b");
+ {
+ test_dim_2(method,
+ 1.2e7, Tensor<1,2>({1,0}),
+ -0.2e-8 );
+ }
+ deallog.pop();
+
+ // Diagonal (equal)
+ deallog.push("Test 2c");
+ {
+ test_dim_2(method,
+ 16.7, Tensor<1,2>({1,0}),
+ 16.7 );
+ }
+ deallog.pop();
+
+ // Non-diagonal
+ deallog.push("Test 2d");
+ {
+ test_dim_2(method,
+ 115.7, Tensor<1,2>({1,1}),
+ 13.6 );
+ }
+ deallog.pop();
+
+ // Non-diagonal (large difference)
+ deallog.push("Test 2e");
+ {
+ const double tol = (method == dealii::ql_implicit_shifts ? 1e-11 : 1e-12);
+ test_dim_2(method,
+ 7.2956e8, Tensor<1,2>({3,2}),
+ -5.284e3, tol );
+ }
+ deallog.pop();
+ }
+
+ // Dim = 3
+ {
+ // Diagonal
+ deallog.push("Test 3a");
+ {
+ test_dim_3(method,
+ 3.6, Tensor<1,3>({1,0,0}),
+ 2.4, Tensor<1,3>({0,1,0}),
+ 1.2);
+ }
+ deallog.pop();
+
+ // Diagonal (large difference)
+ deallog.push("Test 3b");
+ {
+ test_dim_3(method,
+ 1.2e7, Tensor<1,3>({1,0,0}),
+ -0.2e-8, Tensor<1,3>({0,1,0}),
+ -6.5e8);
+ }
+ deallog.pop();
+
+ // Diagonal (2 equal)
+ deallog.push("Test 3c");
+ {
+ test_dim_3(method,
+ 16.7, Tensor<1,3>({1,0,0}),
+ 16.7, Tensor<1,3>({0,1,0}),
+ 1e-6);
+ }
+ deallog.pop();
+
+ // Diagonal (3 equal)
+ deallog.push("Test 3d");
+ {
+ test_dim_3(method,
+ 4.2, Tensor<1,3>({1,0,0}),
+ 4.2, Tensor<1,3>({0,1,0}),
+ 4.2);
+ }
+ deallog.pop();
+
+ // Non-diagonal
+ deallog.push("Test 3e");
+ {
+ test_dim_3(method,
+ 115.7, Tensor<1,3>({1,1,1}),
+ 13.6, Tensor<1,3>({-1,1,-1}),
+ -45.2);
+ }
+ deallog.pop();
+
+ // Non-diagonal (1 large difference)
+ deallog.push("Test 3f");
+ {
+ const double tol = (method == dealii::hybrid ? 1e-9 : (method == dealii::ql_implicit_shifts ? 1e-10 : 5e-11));
+ test_dim_3(method,
+ 7.2956e8, Tensor<1,3>({3,2,5}),
+ -4.856e3, Tensor<1,3>({-0.2,3,1}),
+ -5.284e3, tol);
+ }
+ deallog.pop();
+
+ // Non-diagonal (2 large difference)
+ deallog.push("Test 3g");
+ {
+ const double tol = (method == dealii::hybrid ? 1e-8 : (method == dealii::ql_implicit_shifts ? 1e-7 : 2.5e-10));
+ test_dim_3(method,
+ 9.274e7, Tensor<1,3>({2,-0.7,1.4}),
+ 2.59343, Tensor<1,3>({0.5,-0.22,-1.42}),
+ -5.292e8, tol);
+ }
+ deallog.pop();
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+int main()
+{
+ initlog();
+
+ deallog.push("Hybrid");
+ run_tests(dealii::hybrid);
+ deallog.pop();
+
+ deallog.push("QL");
+ run_tests(dealii::ql_implicit_shifts);
+ deallog.pop();
+//
+ deallog.push("Jacobi");
+ run_tests(dealii::jacobi);
+ deallog.pop();
+}
projects / dealii.git / commitdiff