scientific-computing - GSL gsl_blas_sgemm 矩阵乘法对于 n>140 非常不准确

Question

只是为了使用 gsl_blas_sgemm 函数单独测试一个乘以 150x150 的所有矩阵。我得到的错误多达 38 个。为什么它这么不准确？

Answer 1

你确定你说的对吗？我没有问题：

#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;
}

接着

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