Home > K25104 > LinearSystems > LU_pivot.m

LU_pivot

PURPOSE ^

Computes the LU factorization with partial pivoting

SYNOPSIS ^

function [L,U,P]=LU_pivot(A)

DESCRIPTION ^

 Computes the LU factorization with partial pivoting 
   such that P * A = L * U
 Input arguments:
   A, square matrix
 Output arguments:
   L, square matrix of the same dimensions as A containing the lower
     triangular portion of the LU factorization
   U, square matrix of the same dimensions as A containing the upper
     triangular portion of the LU factorization
   P, square matrix of the same dimensions as A, a permutation matrix
     representing row exchanges due to pivoting

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [L,U,P]=LU_pivot(A)
0002 % Computes the LU factorization with partial pivoting
0003 %   such that P * A = L * U
0004 % Input arguments:
0005 %   A, square matrix
0006 % Output arguments:
0007 %   L, square matrix of the same dimensions as A containing the lower
0008 %     triangular portion of the LU factorization
0009 %   U, square matrix of the same dimensions as A containing the upper
0010 %     triangular portion of the LU factorization
0011 %   P, square matrix of the same dimensions as A, a permutation matrix
0012 %     representing row exchanges due to pivoting
0013 
0014 [n,m]=size(A); % finding the size of A
0015 if n~= m
0016     error('input is not a square matrix');
0017 end
0018 L=eye(n);   % initializing L to the identity matrix
0019 P=L;        % initializing P to the identity matrix
0020 U=A;        % initializing U to be A
0021 for k=1:n   % loop calculates one column of L and one row of U at a time
0022     % Note that U holds in its lower portion a modified portion of A
0023     [~, m]=max(abs(U(k:n,k)));  % for the k-th column find the
0024                                     % maximum element (in modulus)on or
0025                                     % below the diagonal, its position
0026                                     % on or below the diagonal is
0027                                     % returned in m
0028     m=m+k-1; % calculate the row index of the maximum element by
0029              % adding the (k-1) elements above the diagonal
0030     % if m = k, do nothing, because the maximum element is already on
0031     % the diagonal
0032     if m~=k
0033         % swap rows m and k in U
0034         temp=U(k,:);    % temp holds the k-th row of U
0035         U(k,:)=U(m,:);  % assign the m-th row of U to the k-th row of U
0036         U(m,:)=temp;    % the m-th row of U becomes temp
0037         % swap rows m and k in P
0038         temp=P(k,:);
0039         P(k,:)=P(m,:);
0040         P(m,:)=temp;
0041         if k >= 2 % swap rows m and k in the portion of L already
0042                   % calculated
0043             temp=L(k,1:k-1);
0044             L(k,1:k-1)=L(m,1:k-1);
0045             L(m,1:k-1)=temp;
0046         end
0047     end
0048     for j=k+1:n
0049         % if U(k,k) = 0, do nothing, because L is already initialised
0050         % to the identity matrix and thus the k-th column is the k-th
0051         % standard basis vector
0052         if abs(U(k,k)) > 1e-12  % not comparing to zero because of
0053                                 % possible rounding errors
0054         L(j,k)=U(j,k)/U(k,k);   % let L be the k-th column of the current
0055                                 % U scaled by the diagonal element
0056         end
0057         U(j,:)=U(j,:)-L(j,k)*U(k,:); % adust U by subtracting the outer
0058                                      % product of of the k-th column of
0059                                      % L and the k-th row of U
0060     end
0061 end
0062

Generated on Mon 18-Jan-2016 10:25:49 by m2html © 2005