


function [x_sol,rel,boost_counter,parker_matrix_size,relir,counter] = malu(A,b,tol,d,NUM)
%%
% Mixed accelerated LU for solving dense linear systems
%
% Based on paper Mixed accelerated techniques for solving dense linear systems,A.S.Papadakis,E.Gallopoulos
%
% By: A.S.Papadakis 
% email: papadakis@ceid.upatras.gr
%
% INPUT:
%   A:        NxN square matrix
%   b:        Nx1 right hand vector
%   tol: 
%             Boosting will be made only if 
%         1.  the position of the current pivot is multiple 
%             of tol (virtual bad pivot).(uncomment line 63 and comment line 65)
%         2.  the current pivot is under the tol
%   d:        depth of butterflies,(d={1,2} is only supported)
%   NUM:      NUM times stewart on NUM bad pivots, then parker
% 
% OUTPUT:
%   x_sol: solution of Ax=b
%   rel: relative error Ax_sol-b
%   boost_counter: the number of boostings
%   parker_matrix_size: size of matrix in which parker is applied
%
%  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.
% 
Asave = A;
n=length(A);  
X = []; Y = []; 
L = eye(n); U=zeros(n);
k=0;
x2=0;
L22=[];
U22=[];
G=[];
H=[];
AA=[];
boost_counter=0;
parker_matrix_size=0;
relir=-1;
counter=-1;
if NUM<=0
   for j=1:n-1 
     L(j+1:n,j) = A(j+1:n,j)/A(j,j);       
     A(j+1:n,j+1:n) = A(j+1:n,j+1:n)-A(j+1:n,j)/A(j,j)*A(j,j+1:n);
   end
else
for j=1:n-1
    %if(mod(n-j+1,tol)==0)
    
    if (abs(A(j,j)) <= tol)
        
        if(boost_counter==NUM)
             k=j;             
             break;
        end
        boost_counter=boost_counter+1;
        
        Y = [Y;[zeros(1,j-1),eye(1,n-j+1)]];
        X = [X,[zeros(j-1,1);eye(n-j+1,1)]];
       
        
        A(j:n,j:n) = A(j:n,j:n)+X(j:n,boost_counter)*Y(boost_counter,j:n);
    end
    L(j+1:n,j) = A(j+1:n,j)/A(j,j);       
    A(j+1:n,j+1:n) = A(j+1:n,j+1:n)-A(j+1:n,j)/A(j,j)*A(j,j+1:n);
end
end

U = triu(A);

if(boost_counter==0)
   y=L\b;
   x_sol=U\y;
   rel=norm(Asave*x_sol-b);
elseif(k==0)
    
   bb=[b,X];

   y=L\bb;
   zz=U\y;

   z=zz(:,1:1); 
   c1=zz(:,2:end);

   theta=c1*((eye(boost_counter)-Y*c1)\(Y*z));

   x_sol=z+theta;

  rel=norm(Asave*x_sol-b);  
else
  A22=A(k:end,k:end);
  U11=U(1:k-1,1:k-1);
  U12=U(1:k-1,k:end);
  L21=L(k:end,1:k-1);

  bb=[b,X];
  y=L\bb;

  y1=y(1:k-1,:);
  y2=y(k:end,:);
  parker_matrix_size=length(A22);
  [z2,rel1,relir1,counter1] = parker_sol_old(A22,y2,2);
  %  [z2,rel1] = parker_sol(A22,y2,2);


  z1=U11\(y1-U12*z2);

  zz=[z1;z2];

  z=zz(:,1:1);
  c1=zz(:,2:end);

  theta=c1*((eye(boost_counter)-Y*c1)\(Y*z));

  x_sol=z+theta;

  rel=norm(Asave*x_sol-b);

  r=b-Asave*x_sol;
  x_f=x_sol;
  counter=0;
  %{
  while(norm(r)>10^(-12))
  	bb=[r,X];
 	y=L\bb;

  	y1=y(1:k-1,:);
  	y2=y(k:end,:);
  	parker_matrix_size=length(A22);
  	[z2,rel] = parker_sol(A22,y2,2);

  	z1=U11\(y1-U12*z2);

 	zz=[z1;z2];

  	z=zz(:,1:1);
  	c1=zz(:,2:end);

  	theta=c1*((eye(boost_counter)-Y*c1)\(Y*z));

  	d=z+theta;
	
	x_f=x_f+d;
	r=b-Asave*x_f;
	counter=counter+1
 
	relir=norm(r)
  end
  %}
end
end

