行列の積を求めるテストコード

matrix_dot.c

#include <stdint.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>

#define MAT_A_ROW    512
#define MAT_A_COL    512
#define MAT_B_ROW    512
#define MAT_B_COL    512


uint32_t *matrix_create(uint32_t row, uint32_t col)
{
    if ((uint64_t)INT32_MAX <= ((uint64_t)row * (uint64_t)col)) {
        return NULL;
    }
    uint32_t *matrix = NULL;
    matrix = malloc(sizeof(uint32_t) * (row * col));
    if (matrix == NULL) {
        return NULL;
    }
    for(uint32_t i = 0; i < (row * col); i++) {
        // matrix[i] = i + 1;
        matrix[i] = 1;
    }
    return matrix;
}

static inline uint32_t mat_idx(uint32_t row, uint32_t col, uint32_t col_length)
{
    return row * col_length + col;
}

uint32_t *matrix_dot_normal(
    int32_t *matA, uint32_t matA_row, uint32_t matA_col,
    int32_t *matB, uint32_t matB_row, uint32_t matB_col
)
{
    if (matA_col != matB_row) {
        // error
        return NULL;
    }

    uint32_t *matC = NULL;
    uint32_t matC_row = matA_row;
    uint32_t matC_col = matB_col;
    matC = malloc(sizeof(uint32_t) * (matC_row * matC_col));
    if (matC == NULL) {
        return NULL;
    }

    for (uint32_t i = 0; i < matC_row; i++) {
        for (uint32_t j = 0; j < matC_col; j++) {
            matC[mat_idx(i, j, matC_col)] = 0;
            for (uint32_t k = 0; k < matA_col; k++) {
                matC[mat_idx(i, j, matC_col)] += matA[mat_idx(i, k, matA_col)] * matB[mat_idx(k, j, matB_col)];
                // printf("%d = %d x %d\n", matA[mat_idx(i, k, matA_col)] * matB[mat_idx(k, j, matB_col)], matA[mat_idx(i, k, matA_col)], matB[mat_idx(k, j, matB_col)]);
                // printf("C[%d][%d] = %d\n", i, j, matC[mat_idx(i, j, matC_col)]);
            }
        }
    }

    return matC;
}

void matrix_print(uint32_t *mat, uint32_t row, uint32_t col)
{
    for (uint32_t i = 0; i < row; i++) {
        for (uint32_t j = 0; j < col; j++) {
            printf("%ld ", mat[mat_idx(i, j, col)]);
        }
        printf("\n");
    }
}

uint64_t matrix_allsum(uint32_t *mat, uint32_t row, uint32_t col)
{
    uint64_t sum = 0;
    for (uint32_t i = 0; i < row; i++) {
        for (uint32_t j = 0; j < col; j++) {
            sum += (uint64_t)mat[mat_idx(i, j, col)];
        }
    }
    return sum;
}
    

int main(void)
{
    uint32_t *matA = matrix_create(MAT_A_ROW, MAT_A_COL);
    uint32_t *matB = matrix_create(MAT_B_ROW, MAT_B_COL);

    if (matA == NULL || matB == NULL) {
        printf("matA/B is NULL\n");
        return -1;
    }

    clock_t normal_start = clock();
    uint32_t *matC_normal = matrix_dot_normal(matA, MAT_A_ROW, MAT_A_COL, matB, MAT_B_ROW, MAT_B_COL);
    clock_t normal_end = clock();

    if (matC_normal == NULL) {
        printf("matC is NULL\n");
        return -1;
    }

    printf("exec time normal: %fsec\n", (double)(normal_end - normal_start) / CLOCKS_PER_SEC);

    uint64_t sum = matrix_allsum(matC_normal, MAT_A_ROW, MAT_B_COL);
    printf("allsum = %ld\n", sum);

    free(matA);
    free(matB);
    free(matC_normal);
    return 0;
}