Re: [Help-gsl] Strange gsl_blas_dsymv() behaviour

[Korusef]

On Tue, Feb 19, 2008 at 12:48 PM, Frank Reininghaus <
address@hidden> wrote:

> I don't quite understand why you expect only three different entries in
> the result vector of the multiplication. Your matrix J does *NOT* have
> squares with the value 801900 around the diagonal as you said, but 9
> diagonals containing this value on both sides of the main diagonal.
> Therefore, it is quite clear that you get lots of different entries in your
> result vector.
>

Sorry it took me so long to respond. You were right, the example I posted
was wrong. But this one should be correct reproduction of the error, I've
encountered.

/**
 * Demonstrates bug for gsl_blas_dsymv.
 * Multiplying symetric matrix with squares of values around diagonal, zero
diagonal and the rest of the matrix filled with different values,
 * with vector of one and zeroes with ones indexes belonging to the same
square on diagonal.
 * Expected only 3 different values:
 *
 * 1. highest values will be there twice: When all the seven ones are in the
square around diagonal.
 * 2. lower value will there be seven times: When one of the seven ones is
on the diagonal.
 * 3. lowest values, the rest of the vector: When the index is outside the
range of the square on the diagonal.
 *
 * Instead I got 4 different results, output follows:
 *
sum: 63.355728885790576, vec: 63.355728885790583
sum: 63.355728885790576, vec: 63.355728885790583
sum: 63.355728885790576, vec: 63.355728885790583
sum: 63.355728885790576, vec: 63.355728885790583
sum: 63.355728885790576, vec: 63.355728885790583
wrong number of different elements: 4
[-0.663751035816223, 891]
[63.355728885790576, 2]
[63.355728885790583, 5]
[73.915017033422345, 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 / 9) == (i / 9))
      {
        gsl_matrix_set(J, i, j, static_cast<double>(801900 + 891)/ 76027.0);
      }
      else
      {
        gsl_matrix_set(J, i, j, static_cast<double>(-8100 + 891)/ 76027.0);
      }
    }
  }
  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);
  std::cerr << std::fixed;
  std::cerr.precision(15);
  for (Index i = 0; i < 900; ++i)
  {
    double sum = 0;
    for (Index j = 0; j < 900; ++j)
    {
      if (i < j)
      {
        sum += gsl_matrix_get(J, i, j) * gsl_vector_get(X, j);
      }
      else if (i > j)
      {
        sum += gsl_matrix_get(J, j, i) * gsl_vector_get(X, j);
      }
    }
    if (gsl_vector_get(V, i) != sum) std::cerr << "sum: " << sum << ", vec:
" << gsl_vector_get(V, i) << "\n";
  }
  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";
  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;
}

-- 
Zdravi Korusef [Libor Dener]
           (: CauCau :)