// ---------------------------------------------------------------------------- // Numerical diagonalization of 3x3 matrcies // 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 // ---------------------------------------------------------------------------- /* * Changes for IMOD: combined include, removed "inline", declared loop variables * and dimensioned u, q to 3 for compilers on Windows * \$Id\$ */ #include #include #include "dsyev3.h" // Macros #define SQR(x) ((x)*(x)) // x^2 // ---------------------------------------------------------------------------- void dsytrd3(double A[3][3], double Q[3][3], double d[3], double e[2]) // ---------------------------------------------------------------------------- // Reduces a symmetric 3x3 matrix to tridiagonal form by applying // (unitary) Householder transformations: // [ d[0] e[0] ] // A = Q . [ e[0] d[1] e[1] ] . Q^T // [ e[1] d[2] ] // The function accesses only the diagonal and upper triangular parts of // A. The access is read-only. // --------------------------------------------------------------------------- { const int n = 3; int i, j; double u[3], q[3]; double omega, f; double K, h, g; // Initialize Q to the identitity matrix #ifndef EVALS_ONLY for (i=0; i < n; i++) { Q[i][i] = 1.0; for (j=0; j < i; j++) Q[i][j] = Q[j][i] = 0.0; } #endif // Bring first row and column to the desired form h = SQR(A[0][1]) + SQR(A[0][2]); if (A[0][1] > 0) g = -sqrt(h); else g = sqrt(h); e[0] = g; f = g * A[0][1]; u[1] = A[0][1] - g; u[2] = A[0][2]; omega = h - f; if (omega > 0.0) { omega = 1.0 / omega; K = 0.0; for (i=1; i < n; i++) { f = A[1][i] * u[1] + A[i][2] * u[2]; q[i] = omega * f; // p K += u[i] * f; // u* A u } K *= 0.5 * SQR(omega); for (i=1; i < n; i++) q[i] = q[i] - K * u[i]; d[0] = A[0][0]; d[1] = A[1][1] - 2.0*q[1]*u[1]; d[2] = A[2][2] - 2.0*q[2]*u[2]; // Store inverse Householder transformation in Q #ifndef EVALS_ONLY for (j=1; j < n; j++) { f = omega * u[j]; for (i=1; i < n; i++) Q[i][j] = Q[i][j] - f*u[i]; } #endif // Calculate updated A[1][2] and store it in e[1] e[1] = A[1][2] - q[1]*u[2] - u[1]*q[2]; } else { for (i=0; i < n; i++) d[i] = A[i][i]; e[1] = A[1][2]; } }