From 54712f9d44a6d3304ae924cdb8a412e317375ada Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Thu, 27 Jul 2017 10:50:04 +0200 Subject: [PATCH] Add computation of eigenvectors of rank-2 symmetric tensor. This commit adds functions that compute the eigenvectors (and associated eigenvalues) of a rank-2 symmetric tensor using one of threee methods: 1. Transformation to tridiagonal form and using the QL algorithm with implicit shifting. 2. A hybrid algorithm that preferentially uses an analytical approach and falls back to the QL algorithm if the computation is deemed to be inaccurate. 3. The (robust but expensive) Jacobi algorithm. --- ...7JoachimKoppJean-PaulPelteretEsterComellas | 10 + include/deal.II/base/symmetric_tensor.h | 857 +++++++++++++++++- tests/base/symmetric_tensor_41.cc | 311 +++++++ tests/base/symmetric_tensor_41.output | 4 + 4 files changed, 1176 insertions(+), 6 deletions(-) create mode 100644 doc/news/changes/major/20170727JoachimKoppJean-PaulPelteretEsterComellas create mode 100644 tests/base/symmetric_tensor_41.cc create mode 100644 tests/base/symmetric_tensor_41.output diff --git a/doc/news/changes/major/20170727JoachimKoppJean-PaulPelteretEsterComellas b/doc/news/changes/major/20170727JoachimKoppJean-PaulPelteretEsterComellas new file mode 100644 index 0000000000..ad4e93287d --- /dev/null +++ b/doc/news/changes/major/20170727JoachimKoppJean-PaulPelteretEsterComellas @@ -0,0 +1,10 @@ +New: The eigenvectors of a rank-2 symmetric tensor can now be computed using one +of three approaches through the eigenvectors() function. The three algorithms +that have been implemented are: +1. The QL algorithm with implicit shifting. +2. A hybrid algorithm that preferentially uses an analytical algorithm +and falls back to the QL algorithm if the calculations are deemed +inaccurate. +3. The Jacobi algorithm. +
+(Joachim Kopp, Jean-Paul Pelteret, Ester Comellas, 2017/07/27) diff --git a/include/deal.II/base/symmetric_tensor.h b/include/deal.II/base/symmetric_tensor.h index d45427aba5..f4bc00c5a6 100644 --- a/include/deal.II/base/symmetric_tensor.h +++ b/include/deal.II/base/symmetric_tensor.h @@ -2325,7 +2325,10 @@ Number second_invariant (const SymmetricTensor<2,3,Number> &t) /** - * 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 @@ -2340,20 +2343,22 @@ eigenvalues (const SymmetricTensor<2,1,Number> &T) /** - * 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 * the characteristic polynomial - * $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 * subject to round-off errors * 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 @@ -2398,18 +2403,20 @@ eigenvalues (const SymmetricTensor<2,2,Number> &T) /** - * 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 * the characteristic polynomial - * $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 * subject to round-off errors * 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 @@ -2480,6 +2487,844 @@ eigenvalues (const SymmetricTensor<2,3,Number> &T) +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 + void + tridiagonalize (const dealii::SymmetricTensor<2,dim,Number> &A, + dealii::Tensor<2,dim,Number> &Q, + std::array &d, + std::array &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(); + + // 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 u; + for (int i=1; i < dim; i++) + { + u[i] = A[0][i]; + if (i == 1) + u[i] -= g; + } + + std::array 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 + std::array >,dim> + ql_implicit_shifts (const dealii::SymmetricTensor<2,dim,Number> &A) + { + static_assert(numbers::NumberTraits::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 w; + // The off-diagonal elements of the tridiagonal + std::array 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 e; + for (unsigned int i=0; i >,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 >,dim> eig_vals_vecs; + for (unsigned int e=0; e + std::array >,dim> + jacobi (dealii::SymmetricTensor<2,dim,Number> A) + { + static_assert(numbers::NumberTraits::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()); + + // The diagonal elements of the tridiagonal matrix; + // this will ultimately store the eigenvalues + std::array 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 >,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 >,dim> eig_vals_vecs; + for (unsigned int e=0; e + std::array >,2> + hybrid (const dealii::SymmetricTensor<2,2,Number> &A) + { + static_assert(numbers::NumberTraits::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 w = eigenvalues(A); + + std::array >,dim> eig_vals_vecs; + + Number t, u; // Intermediate storage + t = std::abs(w[0]); + for (unsigned int i=1; i 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::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 + std::array >,3> + hybrid (const dealii::SymmetricTensor<2,3,Number> &A) + { + static_assert(numbers::NumberTraits::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 w = eigenvalues(A); + + t = std::abs(w[0]); + for (unsigned int i=1; i 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::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 >,dim> eig_vals_vecs; + for (unsigned int e=0; e + struct SortEigenValuesVectors + { + typedef std::pair > 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 +std::array >,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 +std::array >,2> +eigenvectors (const SymmetricTensor<2,2,Number> &T, + const enum EigenvectorMethod method = ql_implicit_shifts) +{ + std::array >,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 +std::array >,3> +eigenvectors (const SymmetricTensor<2,3,Number> &T, + const enum EigenvectorMethod method = ql_implicit_shifts) +{ + std::array >,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 diff --git a/tests/base/symmetric_tensor_41.cc b/tests/base/symmetric_tensor_41.cc new file mode 100644 index 0000000000..820e8321d8 --- /dev/null +++ b/tests/base/symmetric_tensor_41.cc @@ -0,0 +1,311 @@ +// --------------------------------------------------------------------- +// +// 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 +#include + +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 +bool +is_unit_vector(const Tensor<1,dim> &v) +{ + return std::abs(v.norm() - 1.0) < 1e-9; +} + +template +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 +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(); +} diff --git a/tests/base/symmetric_tensor_41.output b/tests/base/symmetric_tensor_41.output new file mode 100644 index 0000000000..4d7a130701 --- /dev/null +++ b/tests/base/symmetric_tensor_41.output @@ -0,0 +1,4 @@ + +DEAL:Hybrid::OK +DEAL:QL::OK +DEAL:Jacobi::OK -- 2.39.5