自编的matlab的满秩分解函数

矩阵A是m*n,秩为r的，将其分解为P（m*r，秩为r,即列满秩），和Q（r*n，秩为r,即行满秩），满足A=PQ

function [P,Q]=mypq(A)
B=rref(A);
[m,n]=size(A);
P0(1:m,:)=0;
Q0(:,1:n)=0;
for i=1:m
    flag=1;
   for j=1:n
       if B(i,j)==1    
           for k=1:i-1
               if B(k,j)~=0
                   flag=0;
                   break;
               end
           end
           for k=i+1:m
               if B(k,j)~=0
                   flag=0;
                   break;
               end
           end
           if flag==1
               P0=[P0,A(:,j)];
               Q0=[Q0;B(i,:)];
           end
       end    
   end
end
[m1,n1]=size(P0);
[m2,n2]=size(Q0);
P=P0(:,2:n1);
Q=Q0(2:m2,:);