function Pieber

clear; clf; clc;

%% Anfangswerte
N = 60;
M = N;

%  Transportparameter (d,r,k >= 0)
   % Damit A symmetrisch ist, muss k=0 sein.
   % Das CG-Verfahren funktioniert bei pos.def. Matrizen.
   % Fuer r>0 ist A streng diagonal dominant und
   % somit nach Gerschgorin pos.def.
d = 0.01;
k = 1.2;
r = 0.4;

N1   = N + 1;
N4   = round(N/4);
N1N1 = N1*N1;
h    = 1/N;
tau  = 1/M;

V = [zeros(N4,1); ones(N1 - 2*N4,1); zeros(N4,1)];
V = kron(V, V');
V = reshape(V, N1N1, 1);

%% Iterationsmatrix erzeugen
v43 = zeros(N1N1, 1);   % v43 -> (4, ..., 4, 3, --- 4, ..., 4, 3)
v10 = zeros(N1N1, 1);   % v10 -> (1, ..., 1, 0, --- 1, ..., 1, 0)

for i = 1:N1
    v43(((i-1)*N1 + 1):i*N1) = [4*ones(N1-1,1); 3];
    v10(((i-1)*N1 + 1):i*N1) = [ones(N1-1,1); 0];
end

d0  =  d/h^2*v43 + (k/h + r)*ones(N1N1, 1);
dm1 = -(d/h + k)/h*v10;
d1  = -d/h^2*[0; v10(1:end-1)];
dN1 = -d/h^2*ones(N1N1, 1);

A = speye(N1N1) + tau*...
    spdiags([dN1, dm1, d0, d1, dN1], [-N1,-1,0,1,N1], N1N1, N1N1);

clear d0 d1 dm1 dN1 v43 v10
%% PCG

%if k~=0          % Ax = V  -->  A'Ax = A'V
%    B = A.'*A;   % Das wird in LoesSystem beruecksichtigt
%    save daten B
%end

for n = 0:M
    V = LoesSystem(A, V, k, r);
    surf(reshape(V,N1,N1)); axis([1 N1 1 N1 0 1]);
  %  alpha(0.6); grid off;
    view([1 1 1]); drawnow; pause(0.03);
end

% delete daten.mat

end

function x = LoesSystem(A, b, k, r)

if k~=0
    b = A.'*b;
    %load daten
    %A = B;      % Laut Profiler ist A.'*A 
    A = A.'*A;   % schneller als load daten.
end
    
    M = 1./diag(A);
    % Ich speichere bewusst nur die Diagonale,
    % weil ich gesehen habe, dass M*r sehr viel kostet,
    % auch wenn nur die Diagonale besetzt ist.
    % Der Befehl M.*r ist VIEL schneller.

    x   = b;
    r   = b - A*x;
    rho = r'*r;
    i   = 1;
    while sqrt(rho) > 1e-13*norm(b,2) && i<size(A,1)/2
        Mr   = M.*r;
        tau1 = Mr.'*r;
        if i==1
            beta = 0;
            p = Mr;
        else
            beta = tau1/tau0; 
            p    = Mr + beta*p;
        end
        Ap    = A*p;
        alpha = tau1/(p.'*Ap);
        x     = x + alpha*p;
        r     = r - alpha*Ap;
        rho   = r.'*r;
        tau0  = tau1;    % das naechste mal ist tau0 das tau
        i     = i + 1;   % vom vorherigen schritt
    end


end
