function [x] = LMain_elimination(A,b,prec,n)
%%
% 列主消元法求解线性代数方程组 Ax = b的MATLAB实现
% A为待求解方程组的系数矩阵（要求A为非奇异系数矩阵）
% b为方程组右端常数项（列向量）
% prec为计算过程中的控制常数，当某个选主元或完成消元后的系数小于prec时就停止计算
% n为计算过程中的有效数字长度(用vpa函数保留)
%%
% x为求解的结果
%%
fprintf('列主消元法求解线性方程组 Ax = b');
[col,rol] = size(A); % 求方程组系数矩阵的行数和列数,这里col == rol
M = zeros(col,rol); % 用0初始化M矩阵
x = zeros(col,1); % 用0初始化x向量
%%
% 消元计算:对k = 1,2,…,n-1依次计算
for k = 1:col-1
[maxv,r] = max(abs(A(k:n,k))); % 求第k列的最大值，用maxv保存,r记录的是从第k行起的第r行（最大值所在的行号），也就是第r + k - 1行
if (maxv < prec)
disp('主元过小，停止计算');
return;
end
r = r + k - 1; % 最大值所在的行号
A([r k],:) = A([k r],:); % 交换系数矩阵的第k行和第r行
b([r k]) = b([k r]); % 交换常数项的第k个数和第r个数(列向量)
if (vpa(A(k,k),n) == 0) % 如果出现除数为0的情况就不必再计算，跳出此函数
return;
else
M(k+1:col,k) = vpa(vpa(A(k+1:col,k),n)/vpa(A(k,k),n),n);
A(k+1:col,k+1:col) = vpa(vpa(A(k+1:col,k+1:col),n) - vpa(M(k+1:col,k),n)*vpa(A(k,k+1:col),n),n);
b(k+1:col) = vpa(vpa(b(k+1:col),n)-vpa(M(k+1:col,k),n)*vpa(b(k),n),n);
end
end
%%
% 回代求解
if (vpa(A(col,rol),n)~=0) % 如果除数不为0,就计算最后一个未知数的根
x(col) = vpa(vpa(b(col),n)/vpa(A(col,rol),n),n);
end
for k = col-1:-1:1
if (vpa(A(k,k),n) == 0) % 除数为0就跳出函数
return;
else
% 依次回代求解其余的根
x(k) = vpa(vpa(vpa(b(k),n) - vpa(A(k,k+1:col),n)*vpa(x(k+1:col),n),n)/vpa(A(k,k),n),n);
end
end
end
%%