행렬 곱셈 max_active_levels / set_nested

#include <stdio.h>
#include <omp.h>

#define NDIM 1024

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


void mat_mul_ijk(double c[NDIM][NDIM], double a[NDIM][NDIM], double b[NDIM][NDIM],int level_num);

void init_array(double a[NDIM][NDIM])
{
    int i, j;
    for(i = 0; i < NDIM ; i++)
        for(j = 0; j < NDIM ; j++)
            a[i][j] = rand()%100000 + 1;
}
void 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 i, j;
    srand(time(NULL));

    if (argc < 2)
        return 0;
    
    int level_num = atoi(argv[1]);

    if(level_num < 1 || level_num > 3)
    {
        level_num = 1;
    }

    init_array(a);
    init_array(b);
    init_zero(c);

    double t = omp_get_wtime();
    mat_mul_ijk(c, b, a, level_num);
    t = omp_get_wtime() - t;

    printf("%.3f elapsed\n",t);
    return 0;
}

void mat_mul_ijk(double c[NDIM][NDIM], double a[NDIM][NDIM], double b[NDIM][NDIM], int level_num)
{
    int i, j, k;
    double sum;

    omp_set_nested(1);

    omp_set_max_active_levels(level_num);

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