function [b, iflag] = lu_pp_sl( A, b, ipivt );
% 
% LU_PP_SL  computes the solution of a linear system
% A x = b using the LU-decomposition with partial
% pivoting computed by LU.
% 
% Usage:
%        [b, iflag] = lu_pp_sl( A, b, ipivt )
% 
% input:
%        A:     the LU-decomposition of  A  computed by lu_pp
%
%        ipivt: pivot information for the LU-decomposition of  A  
%               computed by lu_pp.m
%   
%        b:     the right hand side b. b can have more than one column. 
%               In this case a system with multiple right hand sides is solved.
%  
% output:
%        b:     the solution of the linear system if iflag = 0
%
%        iflag: error flag
%               iflag = 0  A is nonsingular, solution x is computed.
%               iflag = 1  dimension of A or of b is not correct
%               iflag > 1  zero diagonal element of U
%                          detected in row iflag+1
%  
%
% Matthias Heinkenschloss
% Jan 23, 2001


iflag = 0;

% get size of A and check dimensions
[m,n]    = size(A);
if ( m ~= n | n ~= size(b,1) )
   iflag = 1;
   return
end

% Solve L y = P b  ( b is overwritten by the solution )
 
for k = 1:n-1

%  compute P_k b 
   l = ipivt(k);
   if( k ~= l ) 
       tmp  = b(k,:);
       b(k,:) = b(l,:);
       b(l,:) = tmp;
   end

%  compute M_k b  
   b(k+1:n,:) = b(k+1:n,:) + A(k+1:n,k) * b(k,:);
end


% Solve U x = y  ( b is overwritten by the solution )
 
for k = n:-1:1
    if ( A(k,k) == 0 )
         iflag = k+1;
         return
    end
    b(k,:) = b(k,:)/A(k,k);
    b(1:k-1,:) = b(1:k-1,:) - A(1:k-1,k) * b(k,:);
end

% end of lu_pp_sl