GSL gsl_blas_sgemm matrix multiply highly inaccurate for n>140

Just to test a multiplied a 150x150 matrix of all ones by itself using the gsl_blas_sgemm function. I'm getting errors of up to 38. Why开发者_运维技巧 is it so inaccurate?

---

Are you sure you're calling it right? I have no problems:

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

int main(int argc, char **argv) {
    gsl_matrix_float *a, *b;
    float err, maxerr;
    int n;

    if (argc != 2) {
        fprintf(stderr,"Usage: %s n\n\tWhere n is size of matrix.\n", argv[0]);
        exit(-1);
    }

    n = atoi(argv[1]);
    if (n<2 || n > 1000) {
        n = 100;
        fprintf(stderr,"%s: invalid n = %s; using n = %d.\n",
                argv[0], argv[1], n );
    }

    a = gsl_matrix_float_alloc(n,n);
    b = gsl_matrix_float_alloc(n,n);

    for (int i=0; i<n; i++)
        for (int j=0; j<n; j++)
            gsl_matrix_float_set(a,i,j,1.);

    gsl_blas_sgemm (CblasNoTrans, CblasNoTrans, 1., a, a, 0., b);

    maxerr = 0.;
    for (int i=0; i<n; i++) {
        for (int j=0; j<n; j++) {
            err = gsl_matrix_float_get(b,i,j)-(float)n;
            if (err > maxerr) maxerr = err;
        }
    }

    printf("Maxerr = %f\n", maxerr);

    gsl_matrix_float_free(a);
    gsl_matrix_float_free(b);
    return 0;
}

and then

$ gcc -o gslblas gslblas.c -lgsl -lgslcblas -std=c99
$ ./gslblas 150
Maxerr = 0.000000