A starting point.
Clear[refR,refC,m,f,g];
refR = {x1, x2, x3};
refC = {y1, y2, y3};
m = Array[Subscript[a, ##] &, {3, 3}];
f = Times[#, refR] &@*Map[Normalize[#, Total] &];
g = Transpose[#] &@Times[#, refC] &@*Map[Normalize[#, Total] &]@
Transpose[#] &;
f@m == Table[(m[[i, j]]/Total[m, {2}][[i]])*refR[[i]], {i, 3}, {j, 3}]
g@m == Table[(m[[j, i]]/Total[m][[i]])*refC[[i]], {j, 3}, {i, 3}]
True.
True.
mat = {{21.2, 29, 19.9}, {9, 48.8, 22.4}, {49.8, 62.2, 7.8}};
refR = {75, 82, 125};
refC = {85, 150, 45};
ComposeList[{f, g, f, g, f, g, f, g, f, g, f, g, f, g}, mat]
%== {mat, mat1, mat2, mat3, mat4, mat5, mat6, mat7, mat8, mat9, mat10,
mat11, mat12, mat13, mat14}
True.
Edit
Since refC = {87, 150, 45}; now the above code work.
mat = {{21.2, 29, 19.9}, {9, 48.8, 22.4}, {49.8, 62.2, 7.8}};
refR = {75, 82, 125};
refC = {87, 150, 45};
result=FixedPoint[g@*f, mat]
Total[result, {2}] == refR
Total[result] == refC
{{24.1287, 32.4684, 18.4029}, {9.81314, 52.342, 19.8448}, {53.0581, 65.1896, 6.75229}}
True.
True