Use after free compiler error when compiling with -noPHM and having two DEC_RC_FREE after eachother
When sor is declared inline, gives
/home/thomas/work/issues/use_after_free/v_cycle_rb.c: In function ‘SACf__MAIN__v_cycle__d_X_X__d_X_X__d’:
/home/thomas/work/issues/use_after_free/v_cycle_rb.c:15639:127: error: pointer used after ‘free’ [-Werror=use-after-free]
15639 | inl_14666__flat_1441, (AKS, (NHD, (NUQ, (INT, (GLO, (NON, (NOT, (NDI, (INT, )))))))))), 1, )
| ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~Otherwise, compiles AND computes the correct answer. When doing the same for Full Multigrd, computes the wrong answer.
use StdIO: all;
use Math: all;
use Array: all;
double f(double x, double y)
{
return -2d * sin(x + y);
}
double u(double x, double y)
{
return sin(x + y);
}
/* Interpolation */
inline
double[m, m2], double[m, m2]
prolongate_rb(double[n, n2] x_red, double[n, n2] x_black)
{
x = with {
([0, 0] <= [i, j] < [n, n] step [2, 2]): x_black[i, j / 2];
([0, 1] <= [i, j] < [n, n] step [2, 2]): x_red [i, j / 2];
([1, 0] <= [i, j] < [n, n] step [2, 2]): x_red [i, (n - j - 1) / 2];
([1, 1] <= [i, j] < [n, n] step [2, 2]): x_black[i, (n - j) / 2];
}: genarray([n, n], 0d);
res_black = with {
([0, 0] <= [i, j] < [2 * n, n] step [2, 1]):
x[i / 2, j];
([1, 0] <= [i, j] < [2 * n, n] step [2, 1]):
0.25 * (x[mod([i / 2 , n - 1 - j], [n, n])] +
x[mod([i / 2 , n - j], [n, n])] +
x[mod([i / 2 + 1, n - 1 - j], [n, n])] +
x[mod([i / 2 + 1, n - j], [n, n])]);
}: genarray([2 * n, n], 0d);
res_red = with {
([0, 0] <= [i, j] < [2 * n, n] step [2, 1]):
0.5 * (x[mod([i / 2, j ], [n, n])] +
x[mod([i / 2, j + 1 ], [n, n])]);
([1, 0] <= [i, j] < [2 * n, n] step [2, 1]):
0.5 * (x[mod([i / 2 , n - 1 - j], [n, n])] +
x[mod([i / 2 + 1, n - 1 - j], [n, n])]);
}: genarray([2 * n, n], 0d);
return (res_red, res_black);
}
/* Full-weighting operator: for P the interpolation matrix and d the
* dimension, restrict is P^t / 2^d. */
inline
double[n2, m2], double[n2, m2]
restrict_rb(double[n, m] x_red, double[n, m] x_black)
{
y = with {
([0, 0] <= [i, j] < [n / 2, n / 2]):
(1d * x_black[2 * i, j] +
0.5 * (x_red[mod([2 * i - 1, m - 1 - j], [n, m])] +
x_red[mod([2 * i , j ], [n, m])] +
x_red[mod([2 * i , j - 1 ], [n, m])] +
x_red[mod([2 * i + 1, m - 1 - j], [n, m])]) +
0.25 * (x_black[mod([2 * i - 1, m - 1 - j], [n, m])] +
x_black[mod([2 * i - 1, m - j], [n, m])] +
x_black[mod([2 * i + 1, m - 1 - j], [n, m])] +
x_black[mod([2 * i + 1, m - j], [n, m])])) / 4d;
}: genarray([n / 2, n / 2], 0d);
y_black = {[i, j] -> y[i, 2 * j]
| [0, 0] <= [i, j] < [n / 2, n / 4] step [2, 1];
[i, j] -> y[i, n / 2 - 2 * j - 1]
| [1, 0] <= [i, j] < [n / 2, n / 4] step [2, 1]};
y_red = {[i, j] -> y[i, 2 * j + 1]
| [0, 0] <= [i, j] < [n / 2, n / 4] step [2, 1];
[i, j] -> y[i, n / 2 - 2 * (j + 1)]
| [1, 0] <= [i, j] < [n / 2, n / 4] step [2, 1]};
return (y_red, y_black);
}
double[n, n2], double[n, n2] to_red_black(double[n, n] x)
| (n2 == n / 2)
{
black = {[i, j] -> x[i, 2 * j]
| [0, 0] <= [i, j] < [n, n / 2] step [2, 1];
[i, j] -> x[i, n - 2 * j - 1]
| [1, 0] <= [i, j] < [n, n / 2] step [2, 1]};
red = {[i, j] -> x[i, 2 * j + 1]
| [0, 0] <= [i, j] < [n, n / 2] step [2, 1];
[i, j] -> x[i, n - 2 * (j + 1)]
| [1, 0] <= [i, j] < [n, n / 2] step [2, 1]};
return (red, black);
}
double[n, n] from_red_black(double[n, n2] red, double[n, n2] black)
| (n2 == n / 2)
{
return with {
([0, 0] <= [i, j] < [n, n] step [2, 2]): black[i, j / 2];
([0, 1] <= [i, j] < [n, n] step [2, 2]): red [i, j / 2];
([1, 0] <= [i, j] < [n, n] step [2, 2]): red [i, (n - j - 1) / 2];
([1, 1] <= [i, j] < [n, n] step [2, 2]): black[i, (n - j) / 2];
}: genarray([n, n], 0d);
}
/**
* Weights are Manhattan distance to the central point, so equivalent to
* weights [w[1], w[0], w[1]] in the more general case.
**/
inline
double[n, n2], double[n, n2]
stencil_rb(double[n, n2] x_red, double[n, n2] x_black, double[2] w)
{
res_black = w[0] * x_black +
w[1] * {[i, j] -> x_red[mod([i - 1, n2 - 1 - j], [n, n2])] +
x_red[mod([i , j ], [n, n2])] +
x_red[mod([i , j - 1 ], [n, n2])] +
x_red[mod([i + 1, n2 - 1 - j], [n, n2])]
| [i, j] < [n, n2]};
res_red = w[0] * x_red +
w[1] * {[i, j] -> x_black[mod([i - 1, n2 - 1 - j], [n, n2])] +
x_black[mod([i , j ], [n, n2])] +
x_black[mod([i , j + 1 ], [n, n2])] +
x_black[mod([i + 1, n2 - 1 - j], [n, n2])]
| [i, j] < [n, n2]};
return (res_red, res_black);
}
inline
double[n, n2], double[n, n2]
stencil_rb(double[n, n2] x_red, double[n, n2] x_black, double[3] w)
{
res_black = w[0] * x_black +
w[1] * {[i, j] -> x_red[mod([i - 1, n2 - 1 - j], [n, n2])] +
x_red[mod([i , j ], [n, n2])] +
x_red[mod([i , j - 1 ], [n, n2])] +
x_red[mod([i + 1, n2 - 1 - j], [n, n2])]
| [i, j] < [n, n2]} +
w[2] * {[i, j] -> x_black[mod([i - 1, n2 - 1 - j], [n, n2])] +
x_black[mod([i - 1, n2 - j], [n, n2])] +
x_black[mod([i + 1, n2 - 1 - j], [n, n2])] +
x_black[mod([i + 1, n2 - j], [n, n2])]
| [i, j] < [n, n2]};
res_red = w[0] * x_red +
w[1] * {[i, j] -> x_black[mod([i - 1, n2 - 1 - j], [n, n2])] +
x_black[mod([i , j ], [n, n2])] +
x_black[mod([i , j + 1 ], [n, n2])] +
x_black[mod([i + 1, n2 - 1 - j], [n, n2])]
| [i, j] < [n, n2]} +
w[2] * {[i, j] -> x_red[mod([i - 1, n2 - 1 - j], [n, n2])] +
x_red[mod([i - 1, n2 - 2 - j], [n, n2])] +
x_red[mod([i + 1, n2 - 1 - j], [n, n2])] +
x_red[mod([i + 1, n2 - 2 - j], [n, n2])]
| [i, j] < [n, n2]};
return (res_red, res_black);
}
double L2_rb(double[d:shp] x, double[d:shp] y)
{
return sqrt(sum(x * x) + sum(y * y));
}
/* Turning this noinline removes the error, and the code computes the correct
result. */
inline
double[n, m], double[n, m] sor(double[n, m] u_red, double[n, m] u_black,
double[n, m] f_red, double[n, m] f_black,
double h, double omega)
| (n == 2 * m && m % 2 == 0)
{
update_red = {[i, j] -> u_red[mod([i - 1, m - 1 - j], [n, m])] +
u_red[mod([i , j ], [n, m])] +
u_red[mod([i , j - 1 ], [n, m])] +
u_red[mod([i + 1, m - 1 - j], [n, m])]
| [i, j] < [n, m]};
u_black = (1d - omega) * u_black +
omega / 4d * (update_red - h * h * f_black);
update_black = {[i, j] -> u_black[mod([i - 1, m - 1 - j], [n, m])] +
u_black[mod([i , j ], [n, m])] +
u_black[mod([i , j + 1 ], [n, m])] +
u_black[mod([i + 1, m - 1 - j], [n, m])]
| [i, j] < [n, m]};
u_red = (1d - omega) * u_red +
omega / 4d * (update_black - h * h * f_red);
return (u_red, u_black);
}
inline
double[n, m], double[n, m] sor_solve(double[n, m] f_red, double[n, m] f_black,
double h, int max_iter)
| (n == 2 * m && m % 2 == 0)
{
u_red = {iv -> 0d | iv < [n, m]};
u_black = {iv -> 0d | iv < [n, m]};
pi = 4d * atan(1d);
omega = 2d / (1d + sin(pi * h));
for (t = 0; t < max_iter; t++) {
u_red, u_black = sor(u_red, u_black, f_red, f_black, h, omega);
}
return (u_red, u_black);
}
specialize double[16384, 8192], double[16384, 8192] v_cycle(double[16384, 8192] F_red, double[16384, 8192] F_black, double h);
specialize double[8192, 4096], double[8192, 4096] v_cycle(double[8192, 4096] F_red, double[8192, 4096] F_black, double h);
double[n, n2] , double[n, n2]
v_cycle(double[n, n2] F_red, double[n, n2] F_black, double h)
{
if (n <= 128) {
pi = 4d * atan(1d);
/* spectral_radius approx 1, so log is not stable */
log_spectral_radius = log1p(-2d * sin(pi * h) / (1d + sin(pi * h)));
/* Larger constant because convergence deteriorates significantly if omega
is not exactly optimal, and we do have roundoff errors. */
iter = toi(200d * log(h) / log_spectral_radius);
U_red, U_black = sor_solve(F_red, F_black, h, iter);
} else {
U_red = {iv -> 0d | iv < [n, n / 2]};
U_black = {iv -> 0d | iv < [n, n / 2]};
U_red, U_black = sor(U_red, U_black, F_red, F_black, h, 2d / 3d);
Au_red, Au_black = stencil_rb(U_red, U_black, [-4d, 1d] / (h * h));
r_red = Au_red - F_red;
r_black = Au_black - F_black;
r2h_red, r2h_black = restrict_rb(r_red, r_black);
error2h_red, error2h_black = v_cycle(r2h_red, r2h_black, 2d * h);
error_red, error_black = prolongate_rb(error2h_red, error2h_black);
U_red = U_red - error_red;
U_black = U_black - error_black;
U_red, U_black = sor(U_red, U_black, F_red, F_black, h, 2d / 3d);
}
return (U_red, U_black);
}
int main()
{
pi = 4d * atan(1d);
a = 0d;
b = 2d * pi;
n = 16384;
h = (b - a) / tod(n);
iter = toi(log(tod(n)));
F = {[i, j] -> f(a + tod(i) * h, a + tod(j) * h)
| [i, j] < [n, n]};
F_red, F_black = to_red_black(F);
U_red, U_black = v_cycle(F_red, F_black, h);
Au_red, Au_black = stencil_rb(U_red, U_black, [-4d, 1d] / (h * h));
r_red = F_red - Au_red;
r_black = F_black - Au_black;
for (t = 1; t < iter; t++) {
corr_red, corr_black = v_cycle(r_red, r_black, h);
U_red = U_red + corr_red;
U_black = U_black + corr_black;
Au_red, Au_black = stencil_rb(U_red, U_black, [-4d, 1d] / (h * h));
r_red = F_red - Au_red;
r_black = F_black - Au_black;
}
Utrue = {[i, j] -> u(a + tod(i) * h, a + tod(j) * h)
| [i, j] < [n, n]};
Utrue_red, Utrue_black = to_red_black(Utrue);
printf("Actual relative error %e\n",
L2_rb(U_red - Utrue_red, U_black - Utrue_black) /
L2_rb(Utrue_red, Utrue_black));
return 0;
}