[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] Not correct norm in Jacobi method for eigen-problem
From: |
Nikolay Simakov |
Subject: |
[Bug-gsl] Not correct norm in Jacobi method for eigen-problem |
Date: |
Wed, 09 Mar 2011 11:57:02 -0500 |
User-agent: |
Mozilla/5.0 (X11; U; Linux i686 (x86_64); en-US; rv:1.9.2.15) Gecko/20110303 Lightning/1.0b2 Thunderbird/3.1.9 |
Hi,
I was messing around eigenvalue problems and I found that the norm for
Jacobi method (eigenvalues calculations) is done for whole matrix, but
should be done only for off-diagonal elements. That's very crucial since
this norm used to move out from Jacobi iteration cycle, and that's the
whole idea to make off-diagonal elements zero.
I knew that Jacobi is not fast enough but if it is a simple method and
since it is already in gsl why not to fix it.
in
eigen/jacobi.c
is:
===============================================
inline static double
norm (gsl_matrix * A)
{
size_t i, j, M = A->size1, N = A->size2;
double sum = 0.0, scale = 0.0, ssq = 1.0;
for (i = 0; i < M; i++)
{
for (j = 0; j < N; j++)
{
double Aij = gsl_matrix_get (A, i, j);
if (Aij != 0.0)
{
double ax = fabs (Aij);
if (scale < ax)
{
ssq = 1.0 + ssq * (scale / ax) * (scale / ax);
scale = ax;
}
else
{
ssq += (ax / scale) * (ax / scale);
}
}
}
}
sum = scale * sqrt (ssq);
return sum;
}
should be (based on Algorithm 8.4.3 - Cyclic Jacobi. Golub & Van Loan,
Matrix Computations)
=============================================
inline static double
norm (gsl_matrix * A)
{
size_t i, j, M = A->size1, N = A->size2;
double sum = 0.0, scale = 0.0, ssq = 1.0;
for (i = 0; i < M; i++)
{
for (j = i+1; j < N; j++)
{
double Aij = gsl_matrix_get (A, i, j);
sum+=Aij*Aij;
}
return 2.0*sum;
}
Best,
Nikolay
- [Bug-gsl] Not correct norm in Jacobi method for eigen-problem,
Nikolay Simakov <=