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 mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/bug-gsl
