% **************************************************************
% gls : Solution of a general linear system
% **************************************************************
% BEM-LAP-MAT project
% Matlab/Freemat codes
% Copyright 2008- Stephen Kirkup
% John Tyndall Nuclear Research Institute
% School of Computing Engineering and Physical Sciences
% University of Central Lancashire - www.uclan.ac.uk 
% Westlakes Campus
% Samuel Lindow Building
% West Lakes Science and Technology Park
% Whitehaven
% Cumbria CA24 3JY
% United Kingdom
% smkirkup@uclan.ac.uk
%
% This open source code can be found at
%  www.boundary-element-method.com/freematlab/gl.m
%
%  Issued under the GNU General Public License 2007, see gpl.txt
%
%  Part of the the author's open source BEM packages. 
%  All codes and manuals can be downloaded from 
%  www.boundary-element-method.com
%
% **************************************************************

% gls returns the solution x,y to a problem of the form
%                                A x = B y +c
% where A and B are n by n real matrices and c is a n-vector 
% under the condition(s)
%                         
%             {\alpha}_i x_i + {\beta}_i y_i = f_i   for i=1..n.
%
% Clearly only one of {\alpha}_i or {\beta}_i can be zero for each i.
%
% The method employed involves forming a linear system of the form Cz=d 
% where the n by n matrix C and the vector d can be determined from A,B
% and the {\alpha}_i and {\beta}_i. A standard LU factorisation 
% solution method is then employed to return a solution. From this the actual
% solutions x,y can be determined.

function [x,y]=gls(A,B,c,n,alpha,beta,F)
   gamma=norm(B,inf)/norm(A,inf);
   test=abs(beta)>abs(gamma*alpha);
   
   for (i=1:n)
      if (test(i))
         Fob=F(i)/beta(i);
         aob=alpha(i)/beta(i);
         for (j=1:n)
           c(j)=c(j)+Fob*B(j,i);
           B(j,i)=-aob*B(j,i);
        end
      else
         Foa=F(i)/alpha(i);
         boa=beta(i)/alpha(i);
         for (j=1:n)
            c(j)=c(j)-Foa*A(j,i);
            A(j,i)=-boa*A(j,i);
         end
      end
   end
   
   A=A-B;
   y=A\c';
   
   for (i=1:n)
      if (test(i))
         x(i)=(F(i)-alpha(i)*y(i))/beta(i);
      else
         x(i)=(F(i)-beta(i)*y(i))/alpha(i);
      end
   end
   
   for (i=1:n)
      if (test(i))
         temp=x(i);
         x(i)=y(i);
         y(i)=temp;
      end
   end
   x=x';
   continue