/* * ---------------------------------------------------------------------------- * Numerical diagonalization of 3x3 matrices * Copyright (C) 2006 Joachim Kopp * ---------------------------------------------------------------------------- * This library is free software; you can 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. * * This library is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * Lesser General Public License for more details. * * You should have received a copy of the GNU Lesser General Public * License along with this library; if not, write to the Free Software * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * ---------------------------------------------------------------------------- * ---------------------------------------------------------------------------- SUBROUTINE DSYEVJ3(A, Q, W) * ---------------------------------------------------------------------------- * Calculates the eigenvalues and normalized eigenvectors of a symmetric 3x3 * matrix A using the Jacobi algorithm. * The upper triangular part of A is destroyed during the calculation, * the diagonal elements are read but not destroyed, and the lower * triangular elements are not referenced at all. * ---------------------------------------------------------------------------- * Parameters: * A: The symmetric input matrix * Q: Storage buffer for eigenvectors * W: Storage buffer for eigenvalues * ---------------------------------------------------------------------------- */ #ifndef SWIFT_DSYEVJ3_H #define SWIFT_DSYEVJ3_H /* Standard includes */ #include /* Local includes */ #include "error.h" #include "inline.h" /** * @brief Numerical diagonalization of 3x3 matrices. * * Taken from Koop (2008), International Journal of Modern Physics C, Volume 19, * Issue 03, pp. 523-548. Modified to only return the eigen values. * * @param A The symmetric input matrix. * @param W (return) Storage buffer for eigenvalues */ INLINE static void dsyevj3(double A[3][3], double W[3]) { const int N = 3; /* Initialize W to diag(A) */ for (int X = 0; X < N; X++) W[X] = A[X][X]; /* Calculate SQR(tr(A)) */ double SD = 0.0; for (int X = 0; X < N; X++) SD += fabs(W[X]); SD *= SD; /* Main iteration loop */ for (int I = 0; I < 50; I++) { /* Test for convergence */ double SO = 0.0; for (int X = 0; X < N; X++) for (int Y = X + 1; Y < N; Y++) SO += fabs(A[X][Y]); if (SO == 0.0) return; const double THRESH = (I < 4) ? 0.2 * SO / (N * N) : 0.0; /* Do sweep */ for (int X = 0; X < N; X++) { for (int Y = X + 1; Y < N; Y++) { const double G = 100.0 * (fabs(A[X][Y])); if ((I > 4) && (fabs(W[X]) + G == fabs(W[X])) && (fabs(W[Y]) + G == fabs(W[Y]))) { A[X][Y] = 0.0; } else if (fabs(A[X][Y]) > THRESH) { /* Calculate Jacobi transformation */ const double H = W[Y] - W[X]; double T; if (fabs(H) + G == fabs(H)) T = A[X][Y] / H; else { const double THETA = 0.5 * H / A[X][Y]; if (THETA < 0.0) T = -1.0 / (sqrt(1.0 + THETA * THETA) - THETA); else T = 1.0 / (sqrt(1.0 + THETA * THETA) + THETA); } const double U = 1.0 / sqrt(1.0 + T * T); const double S = T * U; const double Z = T * A[X][Y]; /* Apply Jacobi transformation */ A[X][Y] = 0.0; W[X] = W[X] - Z; W[Y] = W[Y] + Z; for (int R = 0; R < X; R++) { T = A[R][X]; A[R][X] = U * T - S * A[R][Y]; A[R][Y] = S * T + U * A[R][Y]; } for (int R = X + 1; R < Y; R++) { T = A[X][R]; A[X][R] = U * T - S * A[R][Y]; A[R][Y] = S * T + U * A[R][Y]; } for (int R = Y + 1; R < N; R++) { T = A[X][R]; A[X][R] = U * T - S * A[Y][R]; A[Y][R] = S * T + U * A[Y][R]; } /* Update eigenvectors */ /* --- Omitted since only the eigenvalues are desired */ } } } } #ifdef SWIFT_DEBUG_CHECKS /* We should never get here! */ error("No convergence."); #endif } #endif /* SWIFT_DSYEVJ3_H */