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