how to optimize and speed up the multiplication of matrix in c++?

Question

this is optimized implementation of matrix multiplication and this routine performs a matrix multiplication operation. C := C + A * B (where A, B, and C are n-by-n matrices stored in column-major format) On exit, A and B maintain their input values.

void matmul_optimized(int n, int *A, int *B, int *C)
{
    // to the effective bitwise calculation
    // save the matrix as the different type
    int i, j, k;
    int cij;
    for (i = 0; i < n; ++i) {
        for (j = 0; j < n; ++j) {
            cij = C[i + j * n]; // the initialization into C also, add separate additions to the product and sum operations and then record as a separate variable so there is no multiplication
            for (k = 0; k < n; ++k) {
                cij ^= A[i + k * n] & B[k + j * n]; // the multiplication of each terms is expressed by using & operator the addition is done by ^ operator.
            }
            C[i + j * n] = cij; // allocate the final result into C         }
    }
}

how do I more speed up the multiplication of matrix based on above function/method?

this function is tested up to 2048 by 2048 matrix.

the function matmul_optimized is done with matmul.

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

#include "cpucycles.c"
#include "helper_functions.c"
#include "matmul_reference.c"
#include "matmul_optimized.c"

int main()
{
    int i, j;
    int n = 1024; // Number of rows or columns in the square matrices
    int *A, *B;   // Input matrices
    int *C1, *C2; // Output matrices from the reference and optimized implementations

    // Performance and correctness measurement declarations
    long int CLOCK_start, CLOCK_end, CLOCK_total, CLOCK_ref, CLOCK_opt;
    long int COUNTER, REPEAT = 5;
    int difference;
    float speedup;

    // Allocate memory for the matrices
    A = malloc(n * n * sizeof(int));
    B = malloc(n * n * sizeof(int));
    C1 = malloc(n * n * sizeof(int));
    C2 = malloc(n * n * sizeof(int));

    // Fill bits in A, B, C1
    fill(A, n * n);
    fill(B, n * n);
    fill(C1, n * n);

    // Initialize C2 = C1
    for (i = 0; i < n; i++)
        for (j = 0; j < n; j++)
            C2[i * n + j] = C1[i * n + j];

    // Measure performance of the reference implementation
    CLOCK_total = 0;
    for (COUNTER = 0; COUNTER < REPEAT; COUNTER++)
    {
        CLOCK_start = cpucycles();
        matmul_reference(n, A, B, C1);
        CLOCK_end = cpucycles();
        CLOCK_total = CLOCK_total + CLOCK_end - CLOCK_start;
    }
    CLOCK_ref = CLOCK_total / REPEAT;
    printf("n=%d Avg cycle count for reference implementation = %ld\n", n, CLOCK_ref);

    // Measure performance of the optimized implementation
    CLOCK_total = 0;
    for (COUNTER = 0; COUNTER < REPEAT; COUNTER++)
    {
        CLOCK_start = cpucycles();
        matmul_optimized(n, A, B, C2);
        CLOCK_end = cpucycles();
        CLOCK_total = CLOCK_total + CLOCK_end - CLOCK_start;
    }
    CLOCK_opt = CLOCK_total / REPEAT;
    printf("n=%d Avg cycle count for optimized implementation = %ld\n", n, CLOCK_opt);

    speedup = (float)CLOCK_ref / (float)CLOCK_opt;

    // Check correctness by comparing C1 and C2
    difference = 0;
    for (i = 0; i < n; i++)
        for (j = 0; j < n; j++)
            difference = difference + C1[i * n + j] - C2[i * n + j];

    if (difference == 0)
        printf("Speedup factor = %.2f\n", speedup);
    if (difference != 0)
        printf("Reference and optimized implementations do not match\n");

    //print(C2, n);

    free(A);
    free(B);
    free(C1);
    free(C2);
    return 0;
}

Answer 1

You can try algorithm like Strassen or Coppersmith-Winograd and here is also a good example. Or maybe try Parallel computing like future::task or std::thread

Answer 2

Optimizing matrix-matrix multiplication requires careful attention to be paid to a number of issues:

• First, you need to be able to use vector instructions. Only vector instructions can access parallelism inherent in the architecture. So, either your compiler needs to be able to automatically map to vector instructions, or you have to do so by hand, for example by calling the vector intrinsic library for AVX-2 instructions (for x86 architectures).

• Next, you need to pay careful attention to the memory hierarchy. Your performance can easily drop to less than 5% of peak if you don't do this.

• Once you do this right, you will hopefully have broken the computation up into small enough computational chunks that you can also parallelize via OpenMP or pthreads.

A document that carefully steps through what is required can be found at http://www.cs.utexas.edu/users/flame/laff/pfhp/LAFF-On-PfHP.html. (This is very much a work in progress.) At the end of it all, you will have an implementation that gets close to the performance attained by high-performance libraries like Intel's Math Kernel Library (MKL) or the BLAS-like Library Instantiation Software (BLIS).

(And, actually, you CAN then also effectively incorporate Strassen's algorithm. But that is another story, told in Unit 3.5.3 of these notes.)

You may find the following thread relevant: How does BLAS get such extreme performance?