OpenMP. Как распараллелить цикл?

Question

[Екатерина Фролова]

Здравствуйте, помогите, пожалуйста, с распараллеливанием данного цикла на OpenMP. Заранее спасибо.

for (i = 0; i < N; i++)
        {
            result[i] = 0;
            for ( j=RowIndex[i]; j<RowIndex[i + 1]; j++)
            {
                if (i == col[j])
                    result[i] += val[j] * vector[col[j]];
                else
                {
                    result[i] += val[j] * vector[col[j]];
                    result[col[j]] += val[j] * vector[i];
                }
            }
        }

Answer 1

[Roman]

OpenMP GCC ibm.com

Comments

[Екатерина Фролова]

Спасибо, но это я уже читала. Смогла распараллелить, но значения правильные получаю через раз. Поэтому помощи и прошу.

[Roman]

Екатерина Фролова, в коде у вас этого нет, покажите как вы это делаете. Что бы, как минимум, это можно было воспроизвести.

[Екатерина Фролова]

Roman,

#pragma omp parallel for private (i, j) 
        for (i = 0; i < N; i++)
        {
            for ( j=RowIndex[i]; j<RowIndex[i + 1]; j++)
            {
                if (i == col[j])
                    result[i] += val[j] * vector[col[j]];
                else
                {
                    result[i] += val[j] * vector[col[j]];
                    #pragma omp atomic
                    result[col[j]] += val[j] * vector[i];
                }
}

[Roman]

Екатерина Фролова,
а зачем у вас дублируется?
result[i] += val[j] * vector[col[j]];
А в таком варианте

#pragma omp parallel private(i, j)
{
  #pragma omp for 
  for (i = 0; i < N; i++)
  {
    for (j = RowIndex[i]; j < RowIndex[i + 1]; j++)
    {
      #pragma omp atomic 
// lock result[i]
      result[i] += val[j] * vector[col[j]];
// unlock
      if (i != col[j])
      {
        #pragma omp atomic
        result[col[j]] += val[j] * vector[i];
      }
    }
  }
}

[Екатерина Фролова]

Roman, Спасибо, так работает. А возможно как-то избавиться от atomic? На сколько я поняла, его использование сильно замедляет программу.

[Roman]

Екатерина Фролова, разбить на две части
чтобы
result[i] += val[j] * vector[col[j]];
result[col[j]] += val[j] * vector[i];
не пересекались.

[Екатерина Фролова]

Сделала вот так, но теперь вернулась проблема с неверными значениями...

#pragma omp parallel private(i, j)
{
    #pragma omp for 
    for (i = 0; i < N; i++)
    {
        for (j = RowIndex[i]; j < RowIndex[i + 1]; j++)
            result[i] += val[j] * vector[col[j]];    
    }
}

#pragma omp parallel private(i, j)
{
    #pragma omp for 
    for (i = 0; i < N; i++)
    {
        if (i != col[j])
        {
            result[col[j]] += val[j] * vector[i];
        }    
    }
}

[Roman]

Екатерина Фролова, дайте ссылку на ваш код с данными https://pastebin.com

[Екатерина Фролова]

https://pastebin.com/XE4UbeB8
Матрица у меня хранится в файле, в формате mtx (Matrix Market), скачала тут https://sparse.tamu.edu/VDOL/freeFlyingRobot_7

[Roman]

Екатерина Фролова,

code

for(int i = 0; i < N; ++i)
{
  result[i] = 0.0;
}
#pragma omp parallel private(i, j) num_threads(omp_get_max_threads() / 2)
{
    #pragma omp for
		for (i = 0; i < N; i++)
		{
			for (j = RowIndex[i]; j < RowIndex[i + 1]; j++)
			{
        #pragma omp atomic
				result[i] += val[j] * vector[col[j]];
				if (i != col[j])
				{
          #pragma omp atomic
					result[col[j]] += val[j] * vector[i];
				}
			}
		}
}

Проверку нужно делать так

///check multiplication

double eps = 0.00005;
for (int i = 0; i < N; ++i)
{
  if (fabs(result[i] - check_result[i]) < eps)
			continue;
  else
 {
    printf("error-->>%d\n", i);
    printf("check_result_t[%d] = %0.12f\t  result[%d] = %0.12f\n", i, check_result[i], i, result[i]);
    printf("Vectors not match\n");
    break;
 }
}

На сколько я поняла, его использование сильно замедляет программу.

Для вашего случая не сильно.

[Екатерина Фролова]

Благодарю за помощь.