I added OpenACC directives to my red and black Gauss-Seidel solver for the Laplace equation (a simple problem with a hot plate), but GPU-accelerated code is not faster than a processor, even for big problems.
I also wrote a CUDA version, and it is much faster than both (for 512x512, about 2 seconds compared to 25 for CPU and OpenACC).
Can anyone think about the reason for this discrepancy? I understand that CUDA offers the greatest potential speed, but OpenACC should give something better than a processor for big problems (for example, the Jacobi solver for the same problem demonstrated here ).
Here is the relevant code (full working source here ):
#pragma acc data copyin(aP[0:size], aW[0:size], aE[0:size], aS[0:size], aN[0:size], b[0:size]) copy(temp_red[0:size_temp], temp_black[0:size_temp]) // red-black Gauss-Seidel with SOR iteration loop for (iter = 1; iter <= it_max; ++iter) { Real norm_L2 = 0.0; // update red cells #pragma omp parallel for shared(aP, aW, aE, aS, aN, temp_black, temp_red) \ reduction(+:norm_L2) #pragma acc kernels present(aP[0:size], aW[0:size], aE[0:size], aS[0:size], aN[0:size], b[0:size], temp_red[0:size_temp], temp_black[0:size_temp]) #pragma acc loop independent gang vector(4) for (int col = 1; col < NUM + 1; ++col) { #pragma acc loop independent gang vector(64) for (int row = 1; row < (NUM / 2) + 1; ++row) { int ind_red = col * ((NUM / 2) + 2) + row; // local (red) index int ind = 2 * row - (col % 2) - 1 + NUM * (col - 1); // global index #pragma acc cache(aP[ind], b[ind], aW[ind], aE[ind], aS[ind], aN[ind]) Real res = b[ind] + (aW[ind] * temp_black[row + (col - 1) * ((NUM / 2) + 2)] + aE[ind] * temp_black[row + (col + 1) * ((NUM / 2) + 2)] + aS[ind] * temp_black[row - (col % 2) + col * ((NUM / 2) + 2)] + aN[ind] * temp_black[row + ((col + 1) % 2) + col * ((NUM / 2) + 2)]); Real temp_old = temp_red[ind_red]; temp_red[ind_red] = temp_old * (1.0 - omega) + omega * (res / aP[ind]); // calculate residual res = temp_red[ind_red] - temp_old; norm_L2 += (res * res); } // end for row } // end for col // update black cells #pragma omp parallel for shared(aP, aW, aE, aS, aN, temp_black, temp_red) \ reduction(+:norm_L2) #pragma acc kernels present(aP[0:size], aW[0:size], aE[0:size], aS[0:size], aN[0:size], b[0:size], temp_red[0:size_temp], temp_black[0:size_temp]) #pragma acc loop independent gang vector(4) for (int col = 1; col < NUM + 1; ++col) { #pragma acc loop independent gang vector(64) for (int row = 1; row < (NUM / 2) + 1; ++row) { int ind_black = col * ((NUM / 2) + 2) + row; // local (black) index int ind = 2 * row - ((col + 1) % 2) - 1 + NUM * (col - 1); // global index #pragma acc cache(aP[ind], b[ind], aW[ind], aE[ind], aS[ind], aN[ind]) Real res = b[ind] + (aW[ind] * temp_red[row + (col - 1) * ((NUM / 2) + 2)] + aE[ind] * temp_red[row + (col + 1) * ((NUM / 2) + 2)] + aS[ind] * temp_red[row - ((col + 1) % 2) + col * ((NUM / 2) + 2)] + aN[ind] * temp_red[row + (col % 2) + col * ((NUM / 2) + 2)]); Real temp_old = temp_black[ind_black]; temp_black[ind_black] = temp_old * (1.0 - omega) + omega * (res / aP[ind]); // calculate residual res = temp_black[ind_black] - temp_old; norm_L2 += (res * res); } // end for row } // end for col // calculate residual norm_L2 = sqrt(norm_L2 / ((Real)size)); if(iter % 100 == 0) printf("%5d, %0.6f\n", iter, norm_L2); // if tolerance has been reached, end SOR iterations if (norm_L2 < tol) { break; } }