%************************************************************
%        gls.m by www.boundary-element-method.com           *
% ***********************************************************
% Matlab/Freemat/Octave codes
% gls.m Version 1.1 July 2015, tidy up (Version 1, June 2015
% http://www.researchgate.net/profile/Stephen_Kirkup
% School of Computing, Engineering and Physical Sciences
%  University of Central Lancashire
% Issued under the GNU General Public License 2007, see 
%  gpl.txt
%
% Source of the code: 
%  http://www.boundary-element-method.com/mfiles/gls.m
% Source of the user-guide: 
%  http://www.boundary-element-method.com/mfiles/gls_m.pdf
% Test code: 
%  http://www.boundary-element-method.com/mfiles/gls_t.m


function [x, y, xory, L, U, B_gls, P, lfail ] = gls(A, B, c, n, alpha, beta, f)

Tiny = 0.0000000001;
lfail=false;

% Compute the column (1-)norms of the matrices a and b
%  and the 1-norms of the matrices a and b
Anorm = norm(A,1);
Bnorm = norm(B,1);
for j = 1:n
  Acolnorms(j) = norm(A(1:n,j),1);
  Bcolnorms(j) = norm(B(1:n,j),1);
end

% Check that each alpha_i and beta_i are not both zero. If they are then 
%  the subroutine is abandoned
for i = 1:n
 if (abs(alpha(i)) < Tiny * Anorm)
  if (abs(beta(i)) < Tiny * Bnorm)
     'Both alpha(i) and beta(i) are zero. gls is abandoned'
    lfail = true;
    Exit
   end
 end
end

% Set xory. This states whether the x_i or y_i is initially found
for i = 1:n
  if (abs(beta(i)) ^ 2 * Acolnorms(i) > abs(alpha(i)) ^ 2 * Bcolnorms(i)) 
      xory_int(i) = true;
  else
      xory_int(i) = false;
  end
end

A_gls=A;
B_int=B;

% The matrices A and B are then over-written
%  First the columns of the matrices are exchanged when xory_int_j if false
for j = 1:n
   if ~xory_int(j)
       temp = A_gls(1:n, j);
       A_gls(1:n, j) = -B_int(1:n, j);
       B_int(1:n, j) = -temp;
    end
end

%  The B matrix is divided by alpha_j or beta_j, depending on xory_int_j
for j = 1:n
    if xory_int(j)
       alphaorbeta=beta(j);         
    else
       alphaorbeta=alpha(j);        
    end
    B_int(1:n,j)=B_int(1:n,j)/alphaorbeta; 
end

% The terms are subtracted from A
for j = 1:n    
    if xory_int(j)
            A_gls(1:n,j)=A_gls(1:n,j)+alpha(j)*B_int(1:n,j);
    else
            A_gls(1:n,j)=A_gls(1:n,j)+beta(j)*B_int(1:n,j);
    end
end


% Prepare the right hand side for the carrying out forward and backward 
%  substitution
c_int=c+B_int*f;

% Carry out an LU factorisation and solve the system
%  Check for special case that A is diagonal
%   Check to see if A is diagonal
A_diagonal = true;
for i = 1:n
  for j = 1:n
    if ~(i == j) 
        if (abs(A_gls(i, j)) > Tiny * Anorm) 
            A_diagonal = false;
        end
    end
  end
end

%   Solve the system
if (A_diagonal) 
%   Set P to its default
   for i = 1:n
      P(i) = i;
   end
   L=eye(n);
   U=A_gls;
   x_int(i)=c(i)/A_gls(i,i);
   y_int=x_int;
else
%   Carry out an LU factorisation of the matrix A using the lu subroutine.
  [L, U, P]= lu(A_gls);
%   Carry out forward and back substitution
  x_int=LUfbsub(L, U, n, P, c_int);
  y_int=x_int;
end

% Complete the solution
for i = 1:n
  if xory_int(i) 
      y_int(i) = (f(i) - alpha(i) * x_int(i)) / beta(i);
  else
      x_int(i) = (f(i) - beta(i) * y_int(i)) / alpha(i);
  end
end
x=x_int;
y=y_int;
xory=xory_int;
B_gls=B_int;



