w=size(A);
Qm=zeros(w);Rm=zeros(w);
Rm(1,1)=norm(A(:,1));
Qm(:,1)=A(:,1)/Rm(1,1);
for j=2:w
V(:,j)=A(:,j);
for i=1:w-1
Rm(i,j)=Qm(:,i)'*V(:,j);%此处是与也经典算法不同之处
V(:,j)=V(:,j)-Rm(i,j)*Qm(:,i);
end
Rm(j,j)=abs(norm(V(:,j)));
Qm(:,j)=V(:,j)/Rm(j,j);
end
Qm
Rm