[Help-gsl] Strange gsl_blas_dsymv() behaviour

[Korusef]

Hi, I was trying to use GSL for matrix/vector multiplication but strange
error occurred and I received different results than I expected.
I'm not sure whether the mistake is mine or whether it's something inside
GSL.
This is the program that demonstrates the problem:

/**
 * Demonstrates bug for gsl_blas_dsymv.
 * Multiplying symmetric matrix with squares around diagonal of value
801900, zero diagonal and the rest of the matrix filled with -8100
 * with vector of one and zeros with ones indexes belonging to the same
square on diagonal.
 * Multiplying this matrix with vector (720,722,724,725,726,727,728) I
expected only 3 different values:
 *
 * 7 * 801900 = 5613300
 * 6 * 801900 = 4811400
 * 7 * (-8100) = -56700
 *
 * instead of that I got 9 different values
    [-56700.000000, 875]
    [753300.000000, 3]
    [1563300.000000, 3]
    [2373300.000000, 2]
    [3183300.000000, 2]
    [3993300.000000, 3]
    [4803300.000000, 3]
    [4811400.000000, 7]
    [5613300.000000, 2]
 *
 *

g++ (GCC) 4.1.2 (Gentoo 4.1.2)
Copyright (C) 2006 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

g++ gsl_blas_dsymv.cpp -lgsl

*  sci-libs/gsl
      Latest version available: 1.9-r1
      Latest version installed: 1.9-r1
      Size of files: 2,514 kB
      Homepage:      http://www.gnu.org/software/gsl/
      Description:   The GNU Scientific Library
      License:       GPL-2

Linux u-sl 2.6.22-p4smp #4 SMP Mon Feb 11 10:45:43 CET 2008 i686 AMD
Athlon(tm) 64 X2 Dual Core Processor 3800+ AuthenticAMD GNU/Linux

 */
#include <cstdlib>
#include <cassert>
#include <iostream>
#include <iomanip>
#include <map>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

int main(int const, char const* [])
{
  gsl_matrix* J = gsl_matrix_calloc(900, 900);
  typedef int Index;
  for (Index i = 0; i < 900; ++i)
  {
    for (Index j = i + 1; j < 900; ++j)
    {
      if ((j - i) < 9)
      {
        gsl_matrix_set(J, i, j, 801900);
      }
      else
      {
        gsl_matrix_set(J, i, j, -8100);
      }
    }
  }
  gsl_vector* X = gsl_vector_calloc(900);
  //elements around diagonale
  gsl_vector_set(X, 720, 1);
  gsl_vector_set(X, 722, 1);
  gsl_vector_set(X, 724, 1);
  gsl_vector_set(X, 725, 1);
  gsl_vector_set(X, 726, 1);
  gsl_vector_set(X, 727, 1);
  gsl_vector_set(X, 728, 1);
  gsl_vector* V = gsl_vector_calloc(900);
  gsl_blas_dsymv(CblasUpper, 1.0, J, X, 0.0, V);
  typedef std::map<double, Index> TestingMap;
  TestingMap test;
  for (Index i = 0; i < 900; ++i)
  {
    ++(test[gsl_vector_get(V, i)]);
  }
  if (3 != test.size()) std::cerr << "wrong number of different elements: "
<< test.size() << "\n";
  std::cerr << std::fixed;
  for (TestingMap::const_iterator it = test.begin(), itEnd = test.end(); it
!= itEnd; ++it)
  {
    std::cerr << "[" << it->first << ", " << it->second << "]" << "\n";
  }
  gsl_matrix_free(J);
  gsl_vector_free(X);
  gsl_vector_free(V);
  return EXIT_SUCCESS;
}

-- 
Korusef [Libor Dener]