I'm working on LU decomposition of block diagonal matrices using OpenACC.
I get the correct decomposition when I run my code sequentially, whilst when executing it under OpecACC directives I get wrong result when conducting the decomposition.
LU decomposition involves nested loop of that type (see here LUPSolve function):
for (unsigned int i = 0; i < N; i++)
for (unsigned int k = 0; k < i; k++)
It appears that when this type of nested loop is used in a routine seq directive within a parallel region, the device always manages to enter the nested loop even when i=0 (which souldn't be possible because of the k<i condition).
I made a simple code to check it:
#pragma acc routine seq
void test ( int* x, int const n ) {
for (unsigned int i = 0; i < n; i++) {
x[i] = -1;
for (unsigned int k = 0; k < i; k++)
x[i] = k < i;
}
}
int main ( ) {
unsigned const n(4);
unsigned const nb(3);
int x[nb*n];
#pragma acc parallel loop copyout(x[:nb*n])
for (unsigned int b = 0; b < nb; b++)
test(x+b*n,n);
// display x
}
The result I get is this one:
x = 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1,
But the correct one (which I get when I run the code without OpenACC) should be:
x = -1, 1, 1, 1, -1, 1, 1, 1, -1, 1, 1, 1,
I must be doing something wrong because it shouldn't enter the nested loop when i=0...
Plus when I put the loops directly in the parallel region (without using a function call) it does work fine.
for (int i = n-1; i >= 0; i--) { for (int k = i+1; k < n; k++) }And the device doesn't enter the nested loop fori=n-1but if I change this into:for (int i = n-1; i >= 0; i--) { for (int k = 0; k < n-i-1; k++) }The same strange behaviour happens and the device enters the nested loop fori=n-1... It's like0<0is true for nested loop... I just don't know what I'm doing wrong, it must be somewhere in my loop's condition... - Random Tourist