Matrixmultiplikation

/*
 *  Multiplikation von Matrizen
 */

#include <stdio.h>

#define N   300

int main(int argc, char* argv[])
{
  double   a[N][N], b[N][N], c[N][N];
  double   spur, spur1, cth;
  int      i,j,k;
  
  /* Anfangswerte setzen */
  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;
  }
  
  /* Matrixmultiplikation */
  for (i=0; i<N; i++) {
    for (j=0; j<N; j++) {
      c[i][j] = 0.0;
      for (k=0; k<N; k++) {
	c[i][j] += a[i][k] * b[k][j];
      }
    }
  }
  
  /* Testdaten ausgeben */
  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]);
}