function [A_s,B_s,C_s,T,Tinv]=mr_svd(A,B,C,Ts,m)

n=size(A,1);
R = [];
O = [];
for i=0:n-1
  R = [R A^i*B];
  O = [O; C*A^i];
end
[W,S,V]=svd(R*O');
T = sqrt(inv(S))*V'*O;
Tinv = R'*W*sqrt(inv(S));

A_s=T*A*Tinv;
B_s=T*B;
C_s=C*Tinv;

m_n=min(m,size(T,2));

A_s = A_s(1:m_n,1:m_n);
B_s = B_s(1:m_n,:);
C_s = C_s(:,1:m_n);
T = T(1:m_n,:);
Tinv = Tinv(:,1:m_n);

