Fix lu test

[shuds13]

The LU Decomposition test results in Nan's due to the line:

A[i*n+j] /= A[j*n+j];

operating on zeros.

To address the above issue, I have initialized the array based on the original polybench.

• check works with openmp
• check the tulip and other variants.

(Edit: I see https://github.com/cavazos-lab/PolyBench-ACC/blob/master/OpenMP/linear-algebra/solvers/lu/lu.c uses the simplified array. So not sure how that is working. Have not tried running it.)

[shuds13]

OpenMP still got the wrong results.

This OpenMP loop appears incorrect because although the i loop occurs inside the parallel region it is not a parallellized loop and will therefore iterate over all range for every thread.

  #pragma omp parallel default(none) firstprivate(n, A)
  {
  for (int i = 0; i < n; i++) {
    #pragma omp for schedule(static)
    for (int j = 0; j < i; j++) {

This seems to work. But if we use this should probably restructure sequential also to do this. There is not much benefit from the openMP.

for(int k = 0; k < n; k++)
{
    #pragma omp parallel for
    for(int i = k + 1; i < n; i++)
    {
        A[i*n + k] /= A[k*n + k];
        for(int j = k + 1; j < n; j++)
        {
            A[i*n + j] -= A[i*n + k]*A[k*n + j];
        }
    }
}

[shuds13]

An alternative.

static void kernel_lu(int n, double *A)
{
  for (int i = 0; i < n; i++) {

    // 1) Lower part: must be sequential to preserve order
    for (int j = 0; j < i; j++) {
      for (int k = 0; k < j; k++) {
        A[i*n + j] -= A[i*n + k] * A[k*n + j];
      }
      A[i*n + j] /= A[j*n + j];
    }

    // 2) Upper part: can be parallelized
    #pragma omp parallel for
    for (int j = i; j < n; j++) {
      for (int k = 0; k < i; k++) {
        A[i*n + j] -= A[i*n + k] * A[k*n + j];
      }
    }

  }
}