function [x_sol,rel,relir,counter] = parker_sol_old(A,b,d)
%%
% Parker's algorithm for solving dense linera systems
%
% Based on paper Random Butterfly Transformations with Applications in Computational Linear Algebra, D. Stott Parker
%
% By: A.S.Papadakis 
% email: papadakis@ceid.upatras.gr
%
% INPUT:
%   A:        NxN square matrix
%   b:        Nx1 right hand vector
%   d:        depth of butterflies,(d={1,2} is only supported)
% 
% OUTPUT:
%   x_sol: solution of Ax=b
%   rel: relative error Ax_sol-b
%
%  Copyright (C) 2012 A.S.Papadakis
% 
%  This program is free software; you can redistribute it and/or
%  modify it under the terms of the GNU General Public License
%  as published by the Free Software Foundation; either version 2
%  of the License, or (at your option) any later version.
% 
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
% 
n=length(A);
[G,A] = makeButterfly(n,d,A);
H = makeButterfly(length(A),d);

if (G==0)
    x_sol=0;
    rel=inf;
    return;
end

%B=G'*A*H;
B=multi(G,A,H);

L = eye(length(B));
U = B;


for (j = 1:length(B)-1)
    L((j+1):end,j) = U((j+1):end,j)/U(j,j);
    U((j+1):end,j:end) = U((j+1):end,j:end) - L((j+1):end,j)*U(j,j:end);
end


if (length(A)~=n)
        [o,oo]=size(b);
        b=[b;zeros(length(A)-n,oo)];
end

G=comp(G);
H=comp(H);

y=L\(G'*b);


x_hat=U\y;
x=H*x_hat;
   
if (length(A)~=n)
   x=x(1:n,:);
   b=b(1:n,:);
   A=A(1:n,1:n);
end

rel=norm(A*x-b);
x_sol=x;

r=b-A*x_sol;
x_f=x_sol;
counter=0;
relir=0;
%{
while(norm(r)>10^(-13))
	y=L\(G'*r);
	x_hat=U\y;
	d=H*x_hat;

	x_f=x_f+d;
	r=b-A*x_f;
	counter=counter+1;
end

relir=norm(r);
x_sol=x_f;
%}
function B=comp(S)
    if(d==1)
       B11=diag(S(1:length(S)/2));
       B12=diag(S(length(S)/2+1:end));
       B21=B11;
       B22=-B12;
       B=1/sqrt(2)*[B11,B12;B21,B22];
    elseif(d==2)
        d=1;
        B0=comp(S(1:end,2));
        B1=comp(S(1:length(S)/2,1));
        B2=comp(S(length(S)/2+1:end,1));
        W=[B1,zeros(length(B1));zeros(length(B1)),B2];
        B=W*B0;
        d=2;
        
    else
        error('support only for depth=1,2');
    end
        
        
end


function B=multi(G,A,H)
    % input
    % G,H resourcive butterflies of compacted form
    % A matrix
    
    % output
    % B=G'AH
    z=length(A);
    
    if (d==1)
        B=kernel(G,H,A);
    elseif (d==2)
        A11=A(1:z/2,1:z/2);
        A12=A(1:z/2,z/2+1:end);
        A21=A(z/2+1:end,1:z/2);
        A22=A(z/2+1:end,z/2+1:end);
        
        g=G(1:z/2,1);
        h=H(1:z/2,1);
        K11=kernel(g,h,A11);
       
        g=G(1:z/2,1);
        h=H(z/2+1:end,1);
        K12=kernel(g,h,A12);
        
        g=G(z/2+1:end,1);
        h=H(1:z/2,1);
        K21=kernel(g,h,A21);
        
        g=G(z/2+1:end,1);
        h=H(z/2+1:end,1);
        K22=kernel(g,h,A22);
        
        K=[K11,K12;K21,K22];
        
        g=G(1:end,2);
        h=H(1:end,2);
        B=kernel(g,h,K);
    else
        error('Support only for depth=1,2');
    end
end

function B=kernel(g,h,A)
    z=length(A);
    A11=A(1:z/2,1:z/2);
    A12=A(1:z/2,z/2+1:end);
    A21=A(z/2+1:end,1:z/2);
    A22=A(z/2+1:end,z/2+1:end);
    
    C=[A11+A12+A21+A22,A11-A12+A21-A22;A11+A12-A21-A22,A11-A12-A21+A22];
    
    for i=1:z
        for j=1:z
            B(i,j)=1/2*g(i)*C(i,j)*h(j);
        end
    end
 %B=1/2*diag(g)*C*diag(h);

end



function [B,A]=makeButterfly(n,d,A)
    % B resourcive butterfly matrix
    % A augmented matrix with 1s
    % n size of A
    % d depth of B

    if (n<=2^(d-1))
       fprintf('depth of resursive matrix must be less than %d\n',log2(n)+1);
       B=0;
       A=0;
       return;
    end
    if (nargin == 3)
        % if n is not multiple of 2^(d) then we augmenting matrix A with additional 1's on the diagonal
        if(mod(n,2^(d))~=0)
            x= 2^(d) - mod(n,2^(d));
            A11=A;
            A22=diag(ones(x,1));
            A12=zeros(n,x);
            A21=zeros(x,n);

            A=[A11,A12;A21,A22];
            n=length(A);
        end
    end

    B=(1/sqrt(2))*exp((-1/2 + (1/2+1/2).*rand(n,d))/10);
end


end
