HPC - Lecture 8 - Matrix by Matrix multiplication with multithreading (EXERCISE TODO)

HPC - Lecture 8 - Matrix by Matrix multiplication with multithreading

Problem:

• A is a matrix with rows and columns
• B is a matrix with rows and columns
• C is a matrix with rows and columns

It can be observed that each element can be calculated, independently of the others, as the scalar product of the -th row of A and the -th column of B:

• with

In a shared memory MIMD environment (e.g. multicore CPU) it is possible to assign the calculation of each to a thread.

Algorithm 1

Thread Organization - Algorithm 1

Using NT threads, organized as a grid of threads, Algorithm 1 expects each of them to compute, cyclically on rows and columns, elements of the C matrix.

N.B. We assume that and are multiples of NTrow and NTcol.

Algorithm 1

• Algorithm executed by each thread according to the SPMD model.
• Threads are identified by a pair of indices (Idi and Idj).
• The iteration structures with indices ii and jj define the number of elements of matrix C computed by each thread.
• Each element is computed as the dot product (as index varies) of the -th row of A and the -th column of B.

for ii = 0 to N1/Ntrow - 1
    for jj = 0 to N3/Ntcol - 1
        i = Idi + ii * NTrow
        j = Idj + jj * NTcol
        for k = 0 to N2 - 1
            c(i,j) = c(i,j) + a(i,k) * b(k,j)
        endfor
    endfor
endfor

Each thread computes vector-vector dot products.

N.B. For simplicity, we assume and are multiples of Ntrow and NTcol, respectively, so that all threads compute the same number of elements. It can easily be generalized to the case of arbitrary and by assigning a different number of elements to some threads.

Evaluation of algorithm 1

We use Amdahl’s law to evaluate the possible algorithms (P computing units, i.e. the threads).

Let us recall that:

• Speed-Up:
• Absolute Overhead:
• Relative Overhead:

Execution time with threads:

Speed-up derivation: Let us recall that:

Where .

represents the decay factor of the Speed-up.

Total overhead:

Evaluation of algorithm 1 - Performance analysis

For simplicity, let’s consider square matrices of order , and threads.

With thread, you need to access:

• Matrices A and B for reading
• Matrix C for reading and writing.

Remembering that the number of operations is , then:

If , for each element of C, one thread must access a row of A and a column of B for a total of memory accesses and must perform operations.

Since each thread has to compute a total of elements of the matrix C, the execution time with threads is:

The overhead is:

From which we derive:

Even in the absence of critical sections, the parallelization of the problem has introduced an overhead that depends on data access. This is because threads must access the same elements of arrays A and B several times, independently of one another, to calculate their own values of C.

Algorithm 2

A second algorithm can be obtained by observing that each row of the matrix C can be calculated, independently of one another, as the product of the corresponding row of the matrix A and the entire matrix B.

You can assign the calculation of a row to a thread.

Thread Organization - Algorithm 2

Using NT threads arranged along a linear sequence, each of them is then assigned the task of cyclically computing rows of C.

N.B. We assume that is a multiple of NT.

Algorithm 2 (Pseudocode)

• Algorithm executed by each thread according to the SPMD model.
• Threads are identified by an index (ID).
• The iteration structure with index ii defines the rows of matrix C computed by each thread.
• Each row is computed as the dot product (as indices and varies) of the -th row of A and the entire matrix B.

for ii = 0 to N1/NT - 1
    i = ID + ii * NT
    for j = 0 to N3 - 1
        for k = 0 to N2 - 1
            c(i,j) = c(i,j) + a(i,k) * b(k,j)
        endfor
    endfor
endfor

Each thread computes vector-matrix dot products.

Evaluation of algorithm 2

For simplicity, let’s consider square matrices of order , and threads:

With thread, you need to access:

• matrices A and B for reading
• matrix C for reading and writing

Remembering that the number of operations is , then:

If , for each row of C, one thread must access a row of A and the entire matrix B for a total of memory accesses and must perform operations.

Since each thread has to compute a total of elements of the matrix C, the execution time with threads is:

The overhead is:

From which we derive:

This overhead is lower (about half) than that obtained for Algorithm 1.

Important observations

In the two algorithms presented, the overhead depends on the number of memory accesses. In both cases, the number of memory accesses is:

Assuming there is enough cache memory to hold all the problem data, performing matrix multiplication would require a number of memory accesses:

The two studied algorithms introduce an overhead that is considered too high to achieve reasonable speed-up values.

By dividing matrix C into blocks (to be assigned to different threads) and calculating these blocks as the product of matrices performed between blocks of matrix A and blocks of matrix B, we can benefit from a better ratio between the number of memory accesses and the number of operations performed. This leads to a significant reduction in overhead. The hypothesis regarding the size of the cache memory is not restrictive, as it can be overcome through the development of block algorithms.

Algorithm 3

Each block of rows of C can be computed independently of the others as the dot product of the corresponding block of rows of A with the entire matrix B.

We can assign the block calculation of rows to a thread.

Thread organization - Algorithm 3

Using NT threads arranged in order along a linear sequence, each of them is assigned the task of computing one block of rows of C as a matrix-matrix product.

N.B. We assume that is a multiple of NT.

Algorithm 3 Pseudocode:

• Algorithm executed by each thread according to the SPMD model.
• Threads are identified by an index (ID). The index of block is exactly the thread index.
• The iteration structure with index i defines the rows of block computed by the thread with identifier Id.
• This block is computed (as indices vary) as the product of the matrix corresponding to block and the entire matrix B.

for i = ID * N1/NT to (ID + 1) * N1/NT - 1
    for k = 0 to N2 - 1
        for j = 0 to N3 - 1
            c(i,j) = c(i,j) + a(i,k) * b(k,j)
        endfor
    endfor
endfor

Each thread computes 1 matrix-matrix product.

Evaluation of algorithm 3

For square matrices of order , and threads.

With thread, you need to access:

• matrices A and B for reading
• matrix C for reading and writing

Remembering that the number of operations is then:

If , each thread must independently access a row block of A and the entire matrix B in read mode, and a block of C in read and write mode.

Blocks of matrices A and C have rows and columns. Number of operations is . Execution time:

The overhead for algorithm 3 is: from which we derive:

Since usually , this overhead is much lower than that obtained by Algorithm 1 and Algorithm 2. This is achieved thanks to fewer memory accesses relative to the number of operations.

Algorithm 4

Each block of C can be computed independently of the others as the dot product of a row block of A with a column block of B.

We can assign the block calculation to a thread.

Thread organization - Algorithm 4

By arranging the NT threads in a way analogous to the blocks of matrix C, Algorithm 4 requires each of them to compute a block of rows and columns of matrix C.

N.B. We assume that and are integers.

Algorithm 4 (Pseudocode):

• Threads identified by a pair of indices (IDi and IDj). The indices of block is exactly the thread indices.
• Iteration structure with indices i and j defines elements of block computed by the thread with identifier and .
• Block computed (as the indices i,j and k vary) as matrix-matrix product of block and block .

for i = IDi * N1/NTrow to (IDi + 1) * N1/NTrow - 1
    for k = 0 to N2 - 1
        for j = IDj * N3/NTcol to (IDj + 1) * N3/NTcol - 1
            c(i,j) = c(i,j) + a(i,k) * b(k,j)
        endfor
    endfor
endfor

each thread computes 1 matrix-matrix product.

Evaluation of algorithm 3

For simplicity let’s consider square matrices of order , and threads. With thread, you need access to:

• matrices A and B for reading.
• matrix C for reading and writing.

Remembering that the number of operations is , then:

If , each thread must independently access:

• The block , requiring read and write accesses.
• The block , requiring read accesses.
• The block , requiring read accesses.

Execution time:

Overhead:

From which we derive:

The overhead obtained is still lower than that of Algorithm 3. Algorithm 3 represents a special case of Algorithm 4 where and .

---

Part 2: Mat-Mat Product & Multithread: Exercises

Exercise

Implement a multithreaded matrix multiplication version with the following prototype:

void matmatthread (int ldA, int ldB, int ldC, double *A, double *B, double *C, 
                   int N1, int N2, int N3, int db1, int db2, int db3, 
                   int NTROW, int NTCOL) {
    // ...
}

• ldA, ldB, ldC: leading dimensions of the arrays.
• *A, *B, *C: pointers to the matrices.
• N1, N2, N3: matrix dimensions.
• db1, db2, db3: block sizes to be passed to matmatblock.
• NTROW, NTCOL: number of threads in each grid direction.

Important (Call Hierarchy)

The function matmatthread MUST call the function matmatblock.

• matmatthread: Parallel optimization on multicore CPU (via OpenMP).
• matmatblock: Sequential optimization on single core.
• matmat: Basic matrix multiplication.

Exercise

• Compare: matmatblock and matmatthread (Gflops, SpeedUp, Efficiency).
• Matrix size: .
• Block size: .
• Thread configuration:

SpeedUp and Efficiency

Speed-up Results:

• : 1
• : 1.98 (), 1.99 ()
• : 3.92 (), 3.97 ()
• : 7.62 (), 7.73 ()

Efficiency Results:

• : 1
• : 0.99
• : 0.98 (), 0.99 ()
• : 0.95 (), 0.97 ()

Deadline

• Date: Tuesday 05/05.
• Deliverable: Single file matmatthread.c containing:
  • 6 versions of matmat
  • matmatblock
  • matmatthread
• Test with: square matrices , block size=2, .