%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Copyright: Martin Lanser, 2017%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function lu_test_vandermonde(n)

%vandermonde-Matrix
A=fliplr(vander(1:n));

%Speichern, da A ?berschrieben wird
A0=A;

%Vektor f?r implizite Spaltenvertauschung
ind=1:n;

%LU Zerlegung mit impliziter Zeilenpivotsuche
for k=1:n-1
    
    %Pivotsuche
    kmax=k;
    maxx=abs(A(k,ind(k)));
    for i=k+1:n
        if (abs(A(k,ind(i)))>maxx)
            kmax=i;
            maxx=abs(A(k,ind(i)));
        end
    end
    
    %Implizites tauschen
    help=ind(k);
    ind(k)=ind(kmax);
    ind(kmax)=help;
    
    %Elimination mit implizitem Tauschen
    for i=k+1:n
        A(i,ind(k))=A(i,ind(k))/A(k,ind(k));
        for j=k+1:n
            A(i,ind(j))=A(i,ind(j))-A(i,ind(k))*A(k,ind(j));
        end
    end
end

%Das sollte auch das parallele Programm liefern
so_sollte_ergebnis_bei_MPI_aussehen=A

%Um jetzt L und U korrekt aus A zu ziehen, sollten wir einmal das Ergebnis
%explizit permutieren. Also P bauen und A*P rechnen...
P=zeros(n,n);
for i=1:n
    P(ind(i),i)=1;
end

A=A*P;

%L und U rausziehen
L=tril(A);
U=triu(A);

%Eins-Diagonale drauf
for i=1:n
    L(i,i)=1;
end

LxU_ist_vandermonde_permutiert=L*U

test_ob_norm_von_LxU_minus_AxP_gleich_null=norm(L*U-A0*P)