function a=SimpleMatrix(x)%对矩阵实行初等行变换，，化为行最简形矩阵，其特点是：非零行的第一个非零元为1，且这些非零元所在的列的其他元素都为0%x=rand(100);object=x;[m,n]=size(object);for i=1:m temp(1:m,1)=object(1:m,i); if( all(abs(temp(i:m,1))<eps)) continue; end N_zero=find(temp); clear temp; [p,q]=size(N_zero); clear Modulus; if( all(abs(object(i,1:n))<10e-10)) continue; end z=1; for j=1:p T_M=object(N_zero(j,1),i)/object(N_zero(i,1),i); if z==1 Modulus=T_M; z=z+1; continue; end Modulus=[Modulus,T_M]; z=z+1; end %初等行变换：把某一行所有元素的k倍加到另一行对应的元素上去（第j行的k倍加到第i行上，记做ri+k*rj) for j=1:p if j==i continue; end object(N_zero(j,1),1:n)=object(N_zero(j,1),1:n)-object(N_zero(i,1),1:n)*Modulus(j); end L=abs(object)<=10e-10; object(L)=0; %初等行变换：对调两行（对调i,j两行，记做ri<->rj) temp(1:m,1)=object(1:m,i); N_zero=find(temp); clear temp; if N_zero(1)==i object(i,i:n)=object(i,i:n)/object(i,i); continue; end temp(1,1:n)=object(i,1:n); object(i,1:n)=object(N_zero(1,1),1:n); object(N_zero(1,1),1:n)=temp(1,1:n); %初等行变换：以数k～＝0乘某一行中的所有元素（第i行乘k，记作ri*k); object(i,1:n)=object(i,1:n)/object(i,i);endL=abs(object)<=10e-10;object(L)=0;a=object;%特点：n*m矩阵，i行i列或者1或者0，没有处理m-n:m列部分%first step:行排序for i=1:m for j=1:n if abs(a(i,j))>10e-10 sort_seed(i,1)=j; a(i,j:n)=a(i,j:n)/a(i,j); break; end sort_seed(i,1)=n; endendfor i=1:m-1 min=i; for j=i+1:m if sort_seed(j)<sort_seed(min) min=j; temp=sort_seed(i); sort_seed(i)=sort_seed(min); sort_seed(min)=temp; temp(1,1:n)=a(i,1:n); a(i,1:n)=a(min,1:n); a(min,1:n)=temp(1,1:n); end endend%step 2:处理m-n:m列部分N_super=find(sort_seed(:,1)>m);L=abs(a)<10e-10;a(L)=0;if N_super for i=(N_super(1,1)):(m-1) a(i+1,sort_seed(i+1):n)=a(i+1,sort_seed(i+1):n)-a(i,sort_seed(i):n)*a(i+1,sort_seed(i))/a(i,sort_seed(i)); L=abs(a)<=10e-10; a(L)=0; endelse return;end 测试程序：clear;randn('state',0),a=randn(100,120);b=SimpleMatrix(a);c=rref(a);if b==c d=1;else d=[0,0]; e=[0,0]; for i=1:100 for j=1:120 if abs(b(i,j)-c(i,j))>0.0001 d=[d;i,j]; e=[e;b(i,j),c(i,j)]; end end endend 测试表明我做的这个小程序跟matlab提供的rref计算结果相同，但我得程序运算更快点，当然做着个程序前，我还不知道rref，做完后才发现重复工作了。
评论