[Bug-gsl] possible bug in gsl cblas function

[Soumyadip Ghosh]

Hi,

I've been stumbling on some piece of GSL code that has been generating
some spurious results, and to the best of my current understanding, it
happens in a call to a GSL BLAS level 3 function gsl_blas_dsyrk.

I refer to code from GSL 1.6 versions from http://www.gnu.org/software/gsl/. 
>From GSL documentation:

Function: int gsl_blas_dsyrk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t
Trans, double alpha, const gsl_matrix * A, double beta, gsl_matrix *
C)
    These functions compute a rank-k update of the symmetric matrix C,
C = \alpha A A^T + \beta C when Trans is CblasNoTrans and C = \alpha
A^T A + \beta C when Trans is CblasTrans. Since the matrix C is
symmetric only its upper half or lower half need to be stored. When
Uplo is CblasUpper then the upper triangle and diagonal of C are used,
and when Uplo is CblasLower then the lower triangle and diagonal of C
are used.

So, for A an k1xk2 matrix, and C an mxn matrix, we need to check that
 1 -  m==n,  
 2-   n ==k1 if called with CBlasNoTrans, or n == k2  if called with
CBlasTrans.

But the code in ${GSLROOT}/blas/blas.c says 

------------* blas/blas.c code snippet *-------------------
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t K = (Trans == CblasNoTrans) ? A->size2 : A->size1;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
    else if (N != K)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
  
  cblas_dsyrk (CblasRowMajor, Uplo, Trans, INT (N), INT (K), alpha, A->data,
               INT (A->tda), beta, C->data, INT (C->tda));

--------------------end code snippet ----------------
 
So (when CBlasNoTrans is active) we end up with an error if m==n==k2,
which looks contradictory to the stated aim in the documentation.
Assuming the doc is wrong and switched Trans statements, doesn't make
sense either, since then the call cblas_dsyrk produces erroneous
results.

To me the following seems to fix the problem

------------* new blas/blas.c code snippet *-------------------
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t K1 = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  const size_t K2 = (Trans == CblasNoTrans) ? A->size2 : A->size1;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
    else if (N != K1)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
  
  cblas_dsyrk (CblasRowMajor, Uplo, Trans, INT (N), INT (K2), alpha, A->data,
               INT (A->tda), beta, C->data, INT (C->tda));

--------------------end code snippet ----------------

Am I correct in making this change? If not, what may I be missing? If
so, other symmetric rank-k update functions might be affected too.

Thanks!
Soumyadip Ghosh
address@hidden