Matrixmultiplikation mit Blas-Routine

/*
 *  Multiplikation von Matrizen
 */

#include <stdio.h>

#define N   1000


void init(double a[N][N], double b[N][N]) {
  int   i,j;

  for (i=0; i<N; i++) {
    for (j=0; j<N; j++) {
      a[i][j] = i+j;
      b[i][j] = i-j;
    }
    b[i][i] = 1.0;
  }
} 


void test(double c[N][N]) {
  int      i;
  double   spur, spur1, cth;

  spur1 = N * (N-1);
  printf("Spur theoretisch: %lf\n", spur1);

  spur = 0.0;
  for (i=0; i<N; i++) {
    spur += c[i][i];
  }
  printf("Spur            : %lf\n", spur);

  cth = 332334001.0;
  printf("c[0][1] theor.  : %lf\n", cth);
  printf("c[0][1]         : %lf\n", c[0][1]);
}


int main(int argc, char* argv[])
{
  double   a[N][N], b[N][N], c[N][N];
  int      i,j,k;
  
  /* Anfangswerte setzen */
  init(a, b);
  
  /* Matrixmultiplikation */
  {
    int     n     = N;
    char    trans = 'N';
    double  alpha = 1.0;
    double  beta  = 0.0;

    dgemm_(&trans, &trans, &n, &n, &n, &alpha, b, &n, a, &n, &beta, c, &n);
  }

  /* Testdaten ausgeben */
  test(c);
}