|
Packit |
67cb25 |
/* eigen/jacobi.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 2004, 2007 Brian Gough, Gerard Jungman
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_vector.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_matrix.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_eigen.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Algorithm 8.4.3 - Cyclic Jacobi. Golub & Van Loan, Matrix Computations */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static inline double
|
|
Packit |
67cb25 |
symschur2 (gsl_matrix * A, size_t p, size_t q, double *c, double *s)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double Apq = gsl_matrix_get (A, p, q);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (Apq != 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double App = gsl_matrix_get (A, p, p);
|
|
Packit |
67cb25 |
double Aqq = gsl_matrix_get (A, q, q);
|
|
Packit |
67cb25 |
double tau = (Aqq - App) / (2.0 * Apq);
|
|
Packit |
67cb25 |
double t, c1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (tau >= 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
t = 1.0 / (tau + hypot (1.0, tau));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
t = -1.0 / (-tau + hypot (1.0, tau));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
c1 = 1.0 / hypot (1.0, t);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*c = c1;
|
|
Packit |
67cb25 |
*s = t * c1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
*c = 1.0;
|
|
Packit |
67cb25 |
*s = 0.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* reduction in off(A) is 2*(A_pq)^2 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return fabs (Apq);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
inline static void
|
|
Packit |
67cb25 |
apply_jacobi_L (gsl_matrix * A, size_t p, size_t q, double c, double s)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t j;
|
|
Packit |
67cb25 |
const size_t N = A->size2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Apply rotation to matrix A, A' = J^T A */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < N; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double Apj = gsl_matrix_get (A, p, j);
|
|
Packit |
67cb25 |
double Aqj = gsl_matrix_get (A, q, j);
|
|
Packit |
67cb25 |
gsl_matrix_set (A, p, j, Apj * c - Aqj * s);
|
|
Packit |
67cb25 |
gsl_matrix_set (A, q, j, Apj * s + Aqj * c);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
inline static void
|
|
Packit |
67cb25 |
apply_jacobi_R (gsl_matrix * A, size_t p, size_t q, double c, double s)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
const size_t M = A->size1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Apply rotation to matrix A, A' = A J */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < M; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double Aip = gsl_matrix_get (A, i, p);
|
|
Packit |
67cb25 |
double Aiq = gsl_matrix_get (A, i, q);
|
|
Packit |
67cb25 |
gsl_matrix_set (A, i, p, Aip * c - Aiq * s);
|
|
Packit |
67cb25 |
gsl_matrix_set (A, i, q, Aip * s + Aiq * c);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
inline static double
|
|
Packit |
67cb25 |
norm (gsl_matrix * A)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i, j, M = A->size1, N = A->size2;
|
|
Packit |
67cb25 |
double sum = 0.0, scale = 0.0, ssq = 1.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < M; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (j = 0; j < N; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double Aij = gsl_matrix_get (A, i, j);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute norm of off-diagonal elements as per algorithm
|
|
Packit |
67cb25 |
8.4.3 and definition at start of section 8.4.1 */
|
|
Packit |
67cb25 |
if (i == j) continue;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (Aij != 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double ax = fabs (Aij);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (scale < ax)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
ssq = 1.0 + ssq * (scale / ax) * (scale / ax);
|
|
Packit |
67cb25 |
scale = ax;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
ssq += (ax / scale) * (ax / scale);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
sum = scale * sqrt (ssq);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return sum;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_eigen_jacobi (gsl_matrix * a,
|
|
Packit |
67cb25 |
gsl_vector * eval,
|
|
Packit |
67cb25 |
gsl_matrix * evec, unsigned int max_rot, unsigned int *nrot)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i, p, q;
|
|
Packit |
67cb25 |
const size_t M = a->size1, N = a->size2;
|
|
Packit |
67cb25 |
double red, redsum = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (M != N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("eigenproblem requires square matrix", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (M != evec->size1 || M != evec->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("eigenvector matrix must match input matrix", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (M != eval->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("eigenvalue vector must match input matrix", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set_zero (eval);
|
|
Packit |
67cb25 |
gsl_matrix_set_identity (evec);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < max_rot; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double nrm = norm (a);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (nrm == 0.0)
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (p = 0; p < N; p++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (q = p + 1; q < N; q++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double c, s;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
red = symschur2 (a, p, q, &c, &s);
|
|
Packit |
67cb25 |
redsum += red;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute A <- J^T A J */
|
|
Packit |
67cb25 |
apply_jacobi_L (a, p, q, c, s);
|
|
Packit |
67cb25 |
apply_jacobi_R (a, p, q, c, s);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute V <- V J */
|
|
Packit |
67cb25 |
apply_jacobi_R (evec, p, q, c, s);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*nrot = i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (p = 0; p < N; p++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double ep = gsl_matrix_get (a, p, p);
|
|
Packit |
67cb25 |
gsl_vector_set (eval, p, ep);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (i == max_rot)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return GSL_EMAXITER;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_eigen_invert_jacobi (const gsl_matrix * a,
|
|
Packit |
67cb25 |
gsl_matrix * ainv, unsigned int max_rot)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (a->size1 != a->size2 || ainv->size1 != ainv->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("jacobi method requires square matrix", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (a->size1 != ainv->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("inverse matrix must match input matrix", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t n = a->size2;
|
|
Packit |
67cb25 |
size_t i,j,k;
|
|
Packit |
67cb25 |
unsigned int nrot = 0;
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector * eval = gsl_vector_alloc(n);
|
|
Packit |
67cb25 |
gsl_matrix * evec = gsl_matrix_alloc(n, n);
|
|
Packit |
67cb25 |
gsl_matrix * tmp = gsl_matrix_alloc(n, n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_memcpy (tmp, a);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
status = gsl_eigen_jacobi(tmp, eval, evec, max_rot, &nrot);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for(i=0; i
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for(j=0; j
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double ainv_ij = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for(k = 0; k
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double f = 1.0 / gsl_vector_get(eval, k);
|
|
Packit |
67cb25 |
double vik = gsl_matrix_get (evec, i, k);
|
|
Packit |
67cb25 |
double vjk = gsl_matrix_get (evec, j, k);
|
|
Packit |
67cb25 |
ainv_ij += vik * vjk * f;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
gsl_matrix_set (ainv, i, j, ainv_ij);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_free(eval);
|
|
Packit |
67cb25 |
gsl_matrix_free(evec);
|
|
Packit |
67cb25 |
gsl_matrix_free(tmp);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (status)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|