Matlab program for LU Factorization with partial (row) pivoting

	function [L,U,P]=LU_pivot(A)
	% LU factorization with partial (row) pivoting
	% K. Ming Leung, 02/05/03
	[n,n]=size(A);
	L=eye(n); P=L; U=A;
	for k=1:n
	[pivot m]=max(abs(U(k:n,k)));
	m=m+k-1;
	if m~=k
	% interchange rows m and k in U
	temp=U(k,:);
	U(k,:)=U(m,:);
	U(m,:)=temp;
	% interchange rows m and k in P
	temp=P(k,:);
	P(k,:)=P(m,:);
	P(m,:)=temp;
	if k >= 2
	temp=L(k,1:k-1);
	L(k,1:k-1)=L(m,1:k-1);
	L(m,1:k-1)=temp;
	end
	end
	for j=k+1:n
	L(j,k)=U(j,k)/U(k,k);
	U(j,:)=U(j,:)-L(j,k)*U(k,:);
	end
	end

view raw 2013120101.m hosted with ❤ by GitHub

上述代码等价的MATLAB内置函数是[L U P]=lu(A)，此时满足LU=PA。
上述代码若修改成输出P比较好。即对应MATLAB的[L U P]=lu(A,’vector’)。此时满足LU=A(p,:)。

求解AX=B的时候我们需要交换A的行，对应的也需要对B进行相同的行交换。如果写成矩阵乘法的形式的话，LU分解后满足LU=PA，原方程AX=B变为LUX=PB，注意到求解出来的X的行没有发生交换。

万一需要p行变换的逆变换，那么p_back(p)=1:n即可。也就是说p向量比P矩阵更好用。P矩阵意味着需要进行矩阵乘法。

Reference:  http://cis.poly.edu/~mleung/CS3734/s03/ch02/LU_pivot.htm