用“部分”主元高斯消去法进行LU分解，及其Matlab和OpenCV实现

1、介绍
部分选主元中“部分”的意思是我们仅仅在当前列中选择可用的主元。还有一种更完备的选主元方法是所谓的全主元法，它是在未被约化的子矩阵中的所有元素之中选一个最大的，并将它交换到对角元的位置。全主元法虽然在理论上是极为稳定的，但它和部分选主元法比，寻找主元的工作量大大增加，而由于部分选主元法的数值稳定性足以满足实际要求，所以在用高斯消去法求解线性方程组时普遍采用的是部分选主元消去法。
2、伪代码

3、Matlab实现

%% Matlab 2018b——用部分主元高斯消去法进行LU分解：A=LU
A=[1,2,2;4,4,2;4,6,4];
[m,n]=size(A);
I=eye(m,n);     % 单位阵
O=zeros(m,n);   % 0矩阵
L=I;
for i=1:n-1
   P=I;
   M=O;
   [v,p]=max(abs(A(i:m,i)));
   p=p+i-1; % 最大值在A矩阵中的位置（行）
   if i~=p
      P([i,p],:)=P([p,i],:);   % 交换阵 
      A([i,p],:)=A([p,i],:);   % 交换两行
   end
   if A(i,i)==0
       continue;
   else
       for j=i+1:m
          M(j,i)=A(j,i)/A(i,i); 
       end
       for j=i:m
           for k=i+1:n
              A(k,j)=A(k,j)-M(k,i)*A(i,j); 
           end           
       end
   end
   L=L*P'*(I+M);
end
U=A;

4、OpenCV实现

// OpenCV 4.5.0——用部分选主元高斯消去法进行LU分解：A=LU
#include <iostream>
#include <opencv2/opencv.hpp>

using namespace std;
using namespace cv;

int main()
{
	Mat A = (Mat_<float>(3, 3) << 1, 2, 2, 4, 4, 2, 4, 6, 4);
	Size n = A.size();
	Mat I = Mat::eye(n, CV_32FC1);
	Mat O = Mat::zeros(n, CV_32FC1);
	Mat P, L = I.clone();
	Mat M, U = O.clone();
	Mat temp;
	int p;
	for (int i = 0; i < n.width - 1; i++)
	{
		I.copyTo(P);
		O.copyTo(M);
		minMaxIdx(abs(A(Range(i, n.height), Range(i, i + 1))), NULL, NULL, NULL, &p);
		p = p + i;
		if (p!= i)
		{
			// 交换阵
			P.row(p).copyTo(temp);
			P.row(i).copyTo(P.row(p));
			temp.copyTo(P.row(i));
			A = P * A;	// 行交换
		}
		if (A.at<float>(i, i) == 0)
			continue;
		else
		{
			for (int j = i + 1; j < n.height; j++)
				M.at<float>(j, i) = A.at<float>(j, i) / A.at<float>(i, i);
			for (int j = i; j < n.width; j++)
				for (int k = i + 1; k < n.height; k++)
					A.at<float>(k, j) = A.at<float>(k, j) - M.at<float>(k, i)*A.at<float>(i, j);
		}
		L = L * P.t()*(I + M);
	}
	A.copyTo(U);
	cout << "L = " << L << endl;
	cout << "U = " << U << endl;
	getchar();
}

参考文献：《科学计算导论》（第2版）——Michael T. Health 著 张威等 译