LAFF-On Programming for High Performance

Unit 4.3.1 Parallelizing the second loop around the microkernel

Now let's dive right in and parallelize one of the loops. We flipped a coin and it is the second loop around the microkernel that we are going to discuss first:

\lstinputlisting[firstline=73,lastline=79]{github/LAFF-On-PfHPC/Assignments/WeekB/C/GemmFiveLoops.c}

 #pragma omp parallel for
 for ( int j=0; j<n; j+=NR ) {
   int jb = min( NR, n-j );
   LoopOne( m, jb, k, Atilde, &Btilde[ j*k ], &gamma( 0,j ), ldC );
}

lscpu

Architecture:          x86_64
CPU op-mode(s):        32-bit, 64-bit
Byte Order:            Little Endian
CPU(s):                20
On-line CPU(s) list:   0-19
Thread(s) per core:    2
Core(s) per socket:    10
Socket(s):             1
NUMA node(s):          1
Vendor ID:             GenuineIntel
CPU family:            6
Model:                 63
Model name:            Intel(R) Xeon(R) CPU E5-2650 v3 @ 2.30GHz
Stepping:              2
CPU MHz:               1326.273
CPU max MHz:           3000.0000
CPU min MHz:           1200.0000
BogoMIPS:              4594.96
Virtualization:        VT-x
L1d cache:             32K
L1i cache:             32K
L2 cache:              256K
L3 cache:              25600K
NUMA node0 CPU(s):     0-19

sysctl -n hw.physicalcpu

export OMP_NUM_THREADS=4

Unit 4.3.2 Parallelizing the first loop around the microkernel

One can similarly optimize the first loop around the microkernel:

Unit 4.3.3 Parallelizing the third loop around the microkernel

One can similarly optimize the third loop around the microkernel:

Observations:

• Parallelelizing the third loop is a bit trickier. Likely, you started by inserting the #pragma statement and got the wrong answer. The problem is that you are now parallelizing

  
  implemented by

  \begin{center}
  \lstinputlisting[firstline=62,lastline=71]{github/LAFF-On-PfHPC/Assignments/WeekB/C/GemmFiveLoops.c}		
  \end{center}
  

  The problem with a naive parallelization is that all threads pack their block of \(A \) into the same buffer Atilde. There are two solutions to this:

  • Quick and dirty: allocate and free a new buffer in this loop, so that each iteration, and hence each thread, gets a different buffer.

  • Cleaner: In Loop 5 around the microkernel, we currently allocate one buffer for Atilde:

    void LoopThree( int m, int n, int k, double *A, int ldA,
                     double *Btilde, double *C, int ldC )
    {
      double *Atilde =
          ( double * ) malloc( MC * KC * sizeof( double ) );
           
      for ( int i=0; i<m; i+=MC ) {
        int ib = min( MC, m-i );    
        PackBlockA_MCxKC( ib, k, &alpha( i, 0 ), ldA, Atilde );
        LoopTwo( ib, n, k, Atilde, Btilde, &gamma( i,0 ), ldC );
      }
    
      free( Atilde);
    }
    

    What one could do is to allocate a larger buffer here, enough for as many blocks \(\widetilde A\) as there are threads, and then in the third loop around the microkernel one can index into this larger buffer to give each thread its own space. For this, one needs to know about a few OpenMP inquiry routines:

    • omp_get_num_threads() returns the number of theads being used in a parallel section.

    • omp_get_thread_num() returns the thread number of the thread that calls it.

    Now, this is all a bit tricky, because omp_get_num_threads() returns the number of threads being used only in a parallel section. So, to allocate enough space, one would try to change the code in Loop 5 to

    double *Atilde;
    			  
    #pragma omp parallel 
    {
      if ( omp_get_thread_num() == 0 ) 
      malloc( MC * NC * sizeof( double ) * omp_get_num_threads() );
    }
    

    Here the curly brackets indicate the scope of the parallel section. The call to free must then be correspondingly modified so that only thread 0 releases the buffer.

    Now, in Loop 2, one can index into this larger buffer.