1
\$\begingroup\$

Can I improve on this any further. Is there a prettier way of passing the array sizes?

#include <stdio.h>

void matrix_addition(size_t n, int A[n][n], int B[n][n], int C[2*n][2*n], size_t c_start, size_t c_end) {
    for(size_t i = 0; i < n; ++i) {
        for(size_t j = 0; j < n; ++j) {
            C[i+c_start][j+c_end] = A[i][j] + B[i][j];
        }
    }
}

void recursive_matrix_multiply(
        size_t l, size_t n, int A[l][l], int B[l][l], int C[n][n],
        size_t a_r_start, size_t a_r_end, size_t a_c_start, size_t a_c_end, 
        size_t b_r_start, size_t b_r_end, size_t b_c_start, size_t b_c_end
    ) {
    if(n == 1) {
        C[0][0] = A[a_r_start][a_c_start]*B[b_r_start][b_c_start];
    } else {
        int m = n/2;
        int C1[m][m];
        int C2[m][m];

        recursive_matrix_multiply(l, m, A, B, C1, a_r_start, a_r_end-m, a_c_start, a_c_end-m, b_r_start, b_r_end-m, b_c_start, b_c_end-m);
        recursive_matrix_multiply(l, m, A, B, C2, a_r_start, a_r_end-m, a_c_start+m, a_c_end, b_r_start+m, b_r_end, b_c_start, b_c_end-m);
        
        matrix_addition(m, C1, C2, C, 0, 0);
        
        recursive_matrix_multiply(l, m, A, B, C1, a_r_start, a_r_end-m, a_c_start, a_c_end-m, b_r_start, b_r_end-m, b_c_start+m, b_c_end);
        recursive_matrix_multiply(l, m, A, B, C2, a_r_start, a_r_end-m, a_c_start+m, a_c_end, b_r_start+m, b_r_end, b_c_start+m, b_c_end);

        matrix_addition(m, C1, C2, C, 0, m);

        recursive_matrix_multiply(l, m, A, B, C1, a_r_start+m, a_r_end, a_c_start, a_c_end-m, b_r_start, b_r_end-m, b_c_start, b_c_end-m);
        recursive_matrix_multiply(l, m, A, B, C2, a_r_start+m, a_r_end, a_c_start+m, a_c_end, b_r_start+m, b_r_end, b_c_start, b_c_end-m);

        matrix_addition(m, C1, C2, C, m, 0);

        recursive_matrix_multiply(l, m, A, B, C1, a_r_start+m, a_r_end, a_c_start, a_c_end-m, b_r_start, b_r_end-m, b_c_start+m, b_c_end);
        recursive_matrix_multiply(l, m, A, B, C2, a_r_start+m, a_r_end, a_c_start+m, a_c_end, b_r_start+m, b_r_end, b_c_start+m, b_c_end);
    
        matrix_addition(m, C1, C2, C, m, m);
    }
}


int main()
{
    int A[4][4] = {{1, 2, 3, 4}, {4, 3, 2, 1}, { 0, 0, 1, 1 }, {1, 1, 0, 0}};
    int B[4][4] = {{2, 2, 4, 4}, {4, 1, 1, 4}, {1, 0, 1, 0}, {1, 0, 1, 0}};
    int C[4][4];
    
    recursive_matrix_multiply(4, 4, A, B, C, 0, 3, 0, 3, 0, 3, 0, 3);

    for(size_t i = 0; i < 4; ++i) {
        for(size_t j = 0; j < 4; ++j) {
            printf("C[%i][%i] = %i, ", i, j, C[i][j]);
        }
        printf("\n");
    }
    
    return 0;
}

It currently works as intended but I'm looking to beautify this a bit. Even an improvement on variable naming would be appreciated. I'm trying to build up a portfolio so if you saw this from a professional standpoint what would you think?

\$\endgroup\$
2
  • \$\begingroup\$ Is necessary to implement the multiplication in a recursive way ? If not, your matrix multiplication algorithm would take 10 lines. \$\endgroup\$ Commented Jun 19, 2020 at 14:07
  • \$\begingroup\$ @MiguelAvila Yes, I'm avoiding for loops and practicing divide and conquer algorithms. This is really a warm up for Strassens algorithm. \$\endgroup\$ Commented Jun 19, 2020 at 14:15

2 Answers 2

Reset to default
1
\$\begingroup\$

Caution: Assessment of advanced matrix multiplication is non-trivial. See Introduction to Performance Engineering & Matrix Multiplication for starters, in particular the part about cache access patterns and similarities therein between tiling and recursive decomposition on one hand and Strassen&Cie․ on the other.

Public/main items of a program should be self-explanatory, or come with a description.
I'm under the impression that as presented recursive_matrix_multiply() (and consequently matrix_multiply() in the self-answer) only work for square matrices sized the same power of 2.

improvement on variable naming
I can see that telling names aren't short, and only short names keep the length of calls with many parameters and expressions with many primaries from getting out of hand.
I'd much prefer terse names to be introduced with a comment.
I take the names A, B, and C to be conventional for binary operations.
In matrix_addition(), c_start looks plausible, but is it c_ for column or referring to matrix/parameter C?
c_end is worse - I don't see any end it refers to, just the offset for second dimension index in C.
I can live with (I guess I'm still with doxygen)

/// suitable to hold one matrix dimension
typedef size_t md;

void matrix_addition( ///< Let part of C = A + B.
    md n, int A[][n], int B[][n], int C[][2*n],
    md fo, ///< first dimension offset into C
    md so ///< second dimension offset into C
    ) {

I'm avoiding for loops beyond reason -
recursion inside a one nested loop sure is shorter and, in my eyes, more readable.
Taking names from Ivor Denham-Dyson's answer:

    for (md ro = 0 ; ro <= m ; ro += m)
        for (md co = 0 ; co <= m ; co += m) {
            recursive_matrix_multiply(l, m, A, B, C1, ar+ro, ac, br, bc+co);
            recursive_matrix_multiply(l, m, A, B, C2, ar+ro, ac+m, br+m, bc+co);

            matrix_addition_partition(n, m, C1, C2, C, ro, co);
        }
\$\endgroup\$
1
  • \$\begingroup\$ I'm moved on to the suggested text. Such a good resource thanks \$\endgroup\$ Commented Nov 25, 2022 at 19:54
0
\$\begingroup\$

All of the endpoints are redundant. If anyone can add anything else I'd still appreciate it.

#include <stdio.h>

void matrix_addition_partition(size_t n, size_t m, int A[m][m], int B[m][m], int C[n][n], size_t start, size_t end) {
    for(size_t i = 0; i < m; ++i) {
        for(size_t j = 0; j < m; ++j) {
            C[i+start][j+end] = A[i][j] + B[i][j];
        }
    }
}

void recursive_matrix_multiply(
        size_t l, size_t n, int A[l][l], int B[l][l], int C[n][n],
        size_t ar, size_t ac, size_t br, size_t cb
    ) {
    if(n == 1) {
        C[0][0] = A[ar][ac]*B[br][cb];
    } else {
        int m = n>>1;
        int C1[m][m], C2[m][m];

        recursive_matrix_multiply(l, m, A, B, C1, ar, ac, br, cb);
        recursive_matrix_multiply(l, m, A, B, C2, ar, ac+m, br+m, cb);

        matrix_addition_partition(n, m, C1, C2, C, 0, 0);
        
        recursive_matrix_multiply(l, m, A, B, C1, ar, ac, br, cb+m);
        recursive_matrix_multiply(l, m, A, B, C2, ar, ac+m, br+m, cb+m);

        matrix_addition_partition(n, m, C1, C2, C, 0, m);

        recursive_matrix_multiply(l, m, A, B, C1, ar+m, ac, br, cb);
        recursive_matrix_multiply(l, m, A, B, C2, ar+m, ac+m, br+m, cb);

        matrix_addition_partition(n, m, C1, C2, C, m, 0);

        recursive_matrix_multiply(l, m, A, B, C1, ar+m, ac, br, cb+m);
        recursive_matrix_multiply(l, m, A, B, C2, ar+m, ac+m, br+m, cb+m);
    
        matrix_addition_partition(n, m, C1, C2, C, m, m);

    }
}

void matrix_multiply(size_t n, int A[n][n], int B[n][n], int C[n][n]) {
    recursive_matrix_multiply(n, n, A, B, C, 0, 0, 0, 0);
}


int main()
{
    int A[4][4] = {{1, 2, 3, 4}, {4, 3, 2, 1}, {0, 0, 1, 1}, {1, 1, 0, 0}};
    int B[4][4] = {{2, 2, 4, 4}, {4, 1, 1, 4}, {1, 0, 1, 0}, {1, 0, 1, 0}};
    int C[4][4];
    
    matrix_multiply(4, A, B, C);

    for(size_t i = 0; i < 4; ++i) {
        for(size_t j = 0; j < 4; ++j) {
            printf("C[%i][%i] = %i, ", i, j, C[i][j]);
        }
        printf("\n");
    }
    
    return 0;
}
\$\endgroup\$
2
  • \$\begingroup\$ this is hard to parallelize. you can't spawn any threads inside of the recursion as it will create a fork bomb. \$\endgroup\$ Commented Jun 27, 2020 at 14:47
  • \$\begingroup\$ (There was a "finger rot" changing bc to cb.) \$\endgroup\$ Commented Nov 12, 2021 at 8:11

You must log in to answer this question.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.