Hello: I have just noticed that nesting two foralls in a very naive implementation of the over-relation solution of Laplace's equation fails if the compiler's back-end is llvm. In the code below, compiling with chpl --fast laplace-sor-p.chpl and running fails; compiling with chpl --no-llvm --fast laplace-sor-p.chpl works as expected.

I *know* that a better algorithm exists that updates the u's on a "checkers board", so that all new (say "black positions") values are updated using the old (say "white positions") values instead of the completely random update of the code below. I intend to verify if this option passes through llvm without problems.

I am not claiming this is a compiler's bug, since the code below is quite an agressive and thoughtless (naive) implementation, but it may be of interest to know that gcc and llvm backends produce different results for exactly the same source code.

Cheers

Nelson

```
// ===================================================================
// ==> laplace-sor-p: brute force (forall) paralellization of
// over-relaxation with Laplace's equation
// ===================================================================
use Time only stopwatch;
config const epsilon = 1.0e-8; // accuracy
config const omega = 1.95; // relaxation factor
config const Nn = 128; // number of internal points
const dl = 1.0/Nn; // delta l
writef("# Nn = %9i\n", Nn);
writef("# dl = %9.4dr\n", dl);
var u: [0..Nn,0..Nn] real = 0.0; // iterative solution
for i in 0..Nn do { // the initial guess
for j in 0..Nn do {
u[i,j] = IG(i,j);
}
}
var runtime: stopwatch;
var deltam = epsilon; // the norm
var nc = 0; // # of iterations till convergence
runtime.start();
while deltam >= epsilon do { // change to parallel @\label{lin:lapsorp-while}@
deltam = 0.0;
var deltaxm: [1..Nn-1] real ;
forall i in 1..Nn-1 do {
var deltaym: [1..Nn-1] real;
forall j in 1..Nn-1 do {
var uavg = (u[i+1,j]+u[i-1,j]+u[i,j-1]+u[i,j+1])/4.0;
var deltau = omega*(uavg - u[i,j]);
u[i,j] += deltau;
deltau = abs(deltau);
deltaym[j] = deltau;
}
deltaxm[i] = (+ reduce deltaym)/(Nn-1);
}
deltam = (+ reduce deltaxm)/(Nn-1);
writeln("deltam = ", deltam);
nc += 1;
} // done @\label{lin:lapsorp-endwhile}@
runtime.stop();
writeln("nc = ",nc);
var sum = 0.0;
var mad = 0.0;
for i in 0..Nn do {
for j in 0..Nn do {
sum += u[i,j];
var del = abs(u[i,j] - ana(i,j));
mad += del;
}
}
mad /= ((Nn+1)**2);
writeln("uavg = ",sum/((Nn+1)**2));
writeln("mad = ",mad);
writeln("elapsed time = ",runtime.elapsed()," s");
// -------------------------------------------------------------------
// --> initial guess
// -------------------------------------------------------------------
private inline proc IG(i,j: int): real {
const x = i*dl;
const y = j*dl;
if i == 0 then {
return y;
}
else if i == Nn then {
return 1.0;
}
else if j == 0 then {
return x;
}
else if j == Nn then {
return 1.0;
}
else {
return 0.50;
}
}
// -------------------------------------------------------------------
// --> the analytical solution
// -------------------------------------------------------------------
private inline proc ana(i,j): real {
const x = i*dl;
const y = j*dl;
return x + y - x*y;
}
```