[Bug-gsl] Not correct norm in Jacobi method for eigen-problem

[Nikolay Simakov]

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