% Driver for the iterative-solver chapter.
%
% The editable script supplies
%
%     names = method_names()
%     H = run_iter(mid, N, U0, UE, nit, tol)
%
% where H has, in its three rows and one column per iteration, the relative
% residual, the maximum error, and the elapsed milliseconds. Columns past the
% iteration the method stopped at are NaN.
%
% and the driver supplies the starting array U0, which holds the exact
% boundary data on the boundary and zeros inside, and the exact solution UE for
% measurement only. For reference it also solves the same system exactly with
% the sine transform of Chapter 6, which gives both the discretisation error,
% the accuracy no iterative method can beat, and the time a direct solve takes.
p = lp_params();
kind = p.kind;
par = p.par;
N = p.N;
names = method_names();
nm = numel(names);

x = (0:N)' / N;
[X, Y] = meshgrid(x, x);
UE = lp_exact2d(kind, X, Y, par);

U0 = zeros(N+1, N+1);
U0(1, :) = UE(1, :);
U0(N+1, :) = UE(N+1, :);
U0(:, 1) = UE(:, 1);
U0(:, N+1) = UE(:, N+1);

% the exact solution of the discrete system, by the fast Poisson solver
t0 = tic;
Ud = drv_fast(U0, N);
ms_direct = toc(t0) * 1000;
disc_err = max(max(abs(Ud - UE)));

nit = p.nit;
tol = p.tol;
R = NaN(nm, nit+1);
E = NaN(nm, nit+1);
T = NaN(nm, nit+1);
for k = 1:nm
    H = run_iter(k, N, U0, UE, nit, tol);
    R(k, :) = H(1, :);
    E(k, :) = H(2, :);
    T(k, :) = H(3, :);
end

out = struct();
out.names = names;
out.N = N;
out.its = 0:nit;
out.R = lp_rows(R);
out.E = lp_rows(E);
out.T = lp_rows(T);
out.discErr = disc_err;
out.directMs = ms_direct;
out.unknowns = (N-1)^2;
lp_write(out);

function U = drv_fast(U, N)
n = N - 1;
B = zeros(n, n);
B(1, :) = B(1, :) + U(1, 2:N);
B(n, :) = B(n, :) + U(N+1, 2:N);
B(:, 1) = B(:, 1) + U(2:N, 1);
B(:, n) = B(:, n) + U(2:N, N+1);
lam = 2 * (1 - cos(pi * (1:n)' / N));
D = lam * ones(1, n) + ones(n, 1) * lam';
Bh = lp_dst1(lp_dst1(B).').';
U(2:N, 2:N) = lp_dst1(lp_dst1(Bh ./ D).').';
end