/* * 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); }