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
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