Strassen Algorithmus, OMP Version

1. matrix_lib.h

#pragma once

/* Matrix Datenstruktur */
typedef struct  matrix{
    int rows;          /* Anzahl Zeilen/ rows   (m) */
    int cols;          /* Anzahl Spalten / columns (n) */
    double* data;      /* m*n Array in Zeilen Ordnung / row-major */
} matrix_type;

/* Erstelle dynamische Zeilen * Spalten Matrix */
matrix_type* matrix_create(int rows, int cols);

/* Speicher freigeben */
void matrix_free(matrix_type* M);

/* Hilfsfunktionen */
/* Erstellen einer Zufallsmatrix */
void matrix_random(matrix_type* M);

/* Ausgabe Matrix */
void matrix_print(const matrix_type* M, const char* name);

2. matrix_lib.c

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <omp.h>

#include "matrix_lib.h"

matrix_type* matrix_create(int rows, int cols) {
    matrix_type* m = (matrix_type*)malloc(sizeof(matrix_type));
   
    if (!m) {
        fprintf(stderr, "Fehler: malloc matrix_create\n");
        return NULL;
    }
    
    m->rows = rows;
    m->cols = cols;
    m->data = (double*)malloc(rows * cols * sizeof(double));
    if (!m->data) {
        fprintf(stderr, "Fehler: malloc matrix_create.data\n");
        free(m);
        return NULL;
    }
    return m;
}

void matrix_free(matrix_type* m) {
    if (m) {
        if (m->data) {
            free(m->data);
            m->data = NULL;
        }
        free(m);
    }
}

void matrix_random(matrix_type* m) {
    for (int i = 0; i < m->rows * m->cols; i++) {
        m->data[i] = rand() / (double)RAND_MAX;
    }
}

void matrix_print(const matrix_type* m, const char* name) {
    printf("Matrix %s (%dx%d):\n", name, m->rows, m->cols);
    for (int i = 0; i < m->rows; i++) {
        for (int j = 0; j < m->cols; j++) {
            printf("%8.3f ", m->data[i * m->cols + j]);
        }
        putchar('\n');
    }
    putchar('\n');
}

3. strassen.h

#pragma once
#include "matrix_lib.h"

/* Aufruf Strassen Algorithmus */
matrix_type* mul_matrix(const matrix_type* A, const matrix_type* B);

4. strassen.c

#include <stdlib.h>
#include <assert.h>
#include <string.h>
#include <omp.h>
#include "strassen.h"

/*  Deklaration */
matrix_type* omp_multiplikation(const matrix_type* A, const matrix_type* B);

matrix_type* mul_matrix(const matrix_type* A, const matrix_type* B)
{
    assert(A->cols == B->rows);

    int m = A->rows, k = A->cols, n = B->cols;

    /** Groesste Zahl zur Potenz von 2 > max(m, n, k) ermitteln, und in sz speichern */
    int sz = 1;
    while (sz < m || sz < n || sz < k) sz <<= 1;

    matrix_type* Apad = matrix_create(sz, sz);
    matrix_type* Bpad = matrix_create(sz, sz);

    for (int i = 0; i < m; i++) {
        memcpy(Apad->data + i * sz, A->data + i * k, k * sizeof(double));
    }
    for (int i = 0; i < n; i++) {
        memcpy(Bpad->data + i * sz, B->data + i * n, n * sizeof(double));
    }
   
    matrix_type* Cpad = omp_multiplikation(Apad, Bpad);

    matrix_type* C = matrix_create(m, n);
    for (int i = 0; i < m; i++) {
        memcpy(C->data + i * n, Cpad->data + i * sz, n * sizeof(double));
    }

    matrix_free(Apad);
    matrix_free(Bpad);
    matrix_free(Cpad);

    return C;
}

/*  Strassen mit OMP */
matrix_type* omp_multiplikation(const matrix_type* A, const matrix_type* B)
{
    int m = A->rows, k = A->cols, n = B->cols;
    matrix_type* C = matrix_create(m, n);

#pragma omp parallel for schedule(static) collapse(2)
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < n; j++) {
            double sum = 0.0;
#pragma omp simd reduction(+:sum)
            for (int p = 0; p < k; p++) {
                sum += A->data[i * k + p] * B->data[p * n + j];
            }
            C->data[i * n + j] = sum;
        }
    }

    return C;
}

5. main.c

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "matrix_lib.h"
#include "strassen.h"

/*   Programm Strassen - Demo mit Zeitmessung */

int main(void)
{
  
    int n = 512;                      /* Zahl mit der Potenz 2 */

    matrix_type* A = matrix_create(n, n);
    matrix_type* B = matrix_create(n, n);

    matrix_random(A);
    matrix_random(B);

 
#pragma omp parallel
    {
        double start = clock();
        matrix_type* C = mul_matrix(A, B);     /* Strassen Aufruf */
        double end = clock();

        printf("Strassen multiplication von %dx%d Matrizen benoetigt %.3f s\n",n, n, end - start);
        matrix_free(C);
    }

    matrix_free(A);
    matrix_free(B);

    return 0;
}