행렬 곱셈 nest loop 2가지 사용해보기

#include <stdio.h>
#include <omp.h>
#include <stdlib.h>
#include <time.h>
#include "timers.h"

#define NDIM 1024

double a[NDIM][NDIM];
double b[NDIM][NDIM];
double c[NDIM][NDIM];

void matrix_print(double x[NDIM][NDIM])
{
    int i, j;

    for(i = 0; i < NDIM ; i++)
    {
        for(j = 0; j < NDIM ; j++)
            printf("%.2lf ",x[i][j]);
        printf("\n");
    }
}
/*
void matrix_multiply(double a[NDIM][NDIM], double b[NDIM][NDIM], double c[NDIM][NDIM], int num_threads)
{
    int i;
    omp_set_num_threads(num_threads);
#pragma omp parallel
    {
        int i,j,k;
        int tid = omp_get_thread_num();
        int from = (NDIM * tid) / num_threads;
        int to = (NDIM * (tid + 1) ) / num_threads;
        double sum = 0.0;

        for( i = from ; i < to ; i++)
            for(j = 0 ; j < NDIM ; j++)
            {
                sum = 0.0;
                for(k = 0; k < NDIM ; k++)
                {
                    sum += a[i][k] * b[k][j];
                }
                c[i][j] = sum;
            }
    }
}
*/
void matrix_multiply_for_outloop(double a[NDIM][NDIM], double b[NDIM][NDIM], double c[NDIM][NDIM], int num_threads)
{
    int i, j, k;
    double sum = 0.0;
    omp_set_num_threads(num_threads);
#pragma omp parallel 
    {

#pragma omp parallel for private(i, j, k, sum) collapse(2)
        for( i = 0 ; i < NDIM ; i++)
            for(j = 0 ; j < NDIM ; j++)
            {
                sum = 0.0;
                for(k = 0; k < NDIM ; k++)
                {
                    sum += a[i][k] * b[k][j];
                }
                c[i][j] = sum;
            }
    }
}

void matrix_multiply_for_nested_loop(double a[NDIM][NDIM], double b[NDIM][NDIM], double c[NDIM][NDIM], int num_threads)
{
    int i, j, k;
    double sum ;
    omp_set_num_threads(num_threads);
    
    for( i = 0 ; i < NDIM ; i++)
    {

#pragma omp parallel for private(j, k, sum) collapse(2)
            for(j = 0 ; j < NDIM ; j++)
            {
                for(k = 0; k < NDIM ; k++)
                {
                     c[i][j] += a[i][k] * b[k][j];
                }
            }
    }
/*    matrix_print(a);
    printf("====\n");
    matrix_print(b);
    printf("====\n");
    matrix_print(c);
*/
}
void matrix_multiply_for_inner_nested_loop(double a[NDIM][NDIM], double b[NDIM][NDIM], double c[NDIM][NDIM], int num_threads)
{
    int i, j, k;
    double sum = 0.0;
    omp_set_num_threads(num_threads);
#pragma omp parallel 
    {
        int i, j;
        double sum = 0.0;
        for( i = 0 ; i < NDIM ; i++)
            for(j = 0 ; j < NDIM ; j++)
            {
                sum = 0.0;
#pragma omp for private(k)
                for(k = 0; k < NDIM ; k++)
                {
                    sum += a[i][k] * b[k][j];
                }
                c[i][j] = sum;
            }
    }
}


void matrix_init(double a[NDIM][NDIM], double b[NDIM][NDIM])
{
    int i, j;

    srand(time(NULL));
    for(i = 0; i < NDIM ; i++)
        for(j = 0; j < NDIM ; j++)
        {
            a[i][j] = rand()%10 + 1;
            b[i][j] = rand()%10 + 1;
        }
}
void matrix_init_zero(double a[NDIM][NDIM])
{
    int i, j;
    for(i = 0; i < NDIM ; i++)
        for(j = 0; j < NDIM ; j++)
            a[i][j] = 0;
}

int main(int argc, char** argv)
{
    int num_threads = 0;
    
    if(argc > 2)
    {
        num_threads = atoi(argv[1]);
    } else {
        num_threads = 2;
    }

    matrix_init(a, b);

/*    timer_start(1);
    matrix_multiply(a,b,c,num_threads);
    timer_stop(1);

    printf("Time using parallel : %lfsec\n",timer_read(1));
*/    
    timer_start(2);
    matrix_multiply_for_outloop(a,b,c,num_threads);
    timer_stop(2);

    printf("Time using parallel for outloop: %lfsec\n", timer_read(2));

    matrix_init_zero(c);

    timer_start(2);
    matrix_multiply_for_nested_loop(a,b,c,num_threads);
    timer_stop(2);

    printf("Time using parallel for nested loop: %lfsec\n", timer_read(2));

    matrix_init_zero(c);

    timer_start(2);
    matrix_multiply_for_inner_nested_loop(a,b,c,num_threads);
    timer_stop(2);

    printf("Time using parallel for inner nested loop: %lfsec\n", timer_read(2));

    return 0;
}