GSL gsl_blas_sgemm matrice multiplier fortement inexact pour n> 140

Question

Juste pour tester une matrice multipliée par 150x150 de tous les uns par elle-même en utilisant la fonction gsl_blas_sgemm. Je reçois des erreurs allant jusqu'à 38. Pourquoi est-ce si inexact?

Answer 1

Êtes-vous sûr de bien l'appeler? Je n'ai aucun problème:

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

puis

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