% ----------------------------------------------------------------------
%          function LUfbsub.m by www.numerical-methods.com           |
% ----------------------------------------------------------------------
%
% Function LUfbsub returns the solution of a real matrix-vector system
%  of equations
%                              L U x = P b
%
%  with the matrix being a lower-triangular matrix, the matrix U being
%  an upper triangular matrix and P being a permutation matrix. 
%
% The matrics L and U often arise from the LU factorisation of a matrix
%  A.  That is the matrices P, L and U are computed so that
%
%                             PA  = LU
%
%  and in Matlab or Freemat the matrices L,U and P can be determined
%  for a matrix A by using for example the command [L, U, P]=lu(A).
%
% function [x]=LUfbsub(L, U, n, P, b)' Sub LUfbsub(a, n, perm, b)
%  real/complex   L: the nxn lower-triangular matrix 
%  real/complex   U: the nxn upper-triangular matrix 
%  integer        n: the dimension of the matrix/vector
%                 P: an nxn permutation matrix
%  real/complex       b: the vector b
%
% As a result running the function, x contains the solution.
%
% Source of the code: http://www.numerical-methods.com/mfiles/LUfbsub.m
% Source of the user-guide: 
%  http://www.numerical-methods.com/mfiles/LUfbsub_m.htm
%
% Licence: This is 'open source'; the software may be used and applied
%  within other systems as long as its provenance is appropriately
%  acknowledged. See the GNU Licence http://www.gnu.org/licenses/lgpl.txt
%  for more information or contact webmaster@numerical-methods.com

function [x]=LUfbsub(L, U, n, P, b)
    
% Initialise temp as b with rows exchanged as indicated by the
%  the permutation matrix P
temp=P*b;

% Solve Ly=Pb by forward substitution (temp stores y on completion)
% It is presumed that matlab can detect that L is lower triangular
temp=L\temp;

% Solve Ux=y by back substitution (temp stores the solution 'x' on
%  completion)
% It is presumed that matlab can detect that U is upper triangular
temp=U\temp;

x=temp;

