A straightfoward method to check the matrix product you asked for: define your matrix by blocks using ArrayFlatten:

A = Table[RandomReal[], {i, 3}, {j, 3}];
b = List /@ Table[RandomReal[], {i, 3}]; (*So it is a matrix*)
M = ArrayFlatten[{{A, b}}];
ret = ArrayFlatten[{{Inverse[A], -Inverse[A].b}, {0, 1}}]; (*Notice numbers represents matrices for ArrayFlatten*)
M.ret // Chop // MatrixForm

Which returns the 3x4 identity matrix we were expecting.

I guess I should state that Mathematica is doing it right. Truncation error is not avoidable with finite precision arithmetic. As for the problem at hand, there are a couple of things that seem to be amiss.

(1) Create the augmented matrix {A|b] correctly. This can be done with M = Join[A, List /@ b, 2]. Your version transposes A in the process.

(2) Now just use RowReduce on it. No inverses, no multiplications.

Here is an example.

SeedRandom[1111];
aA = Table[RandomReal[], {i, 3}, {j, 3}]
b = Table[RandomReal[], {i, 3}]
mat = Join[aA, List /@ b, 2]
rred = RowReduce[mat]

Out[434]= {{0.0777960697813, 0.556265936296, 
  0.807681686834}, {0.527982841158, 0.555876762387, 
  0.981090055861}, {0.907596678267, 0.364390115148, 0.722450567448}}

Out[435]= {0.420186553434, 0.1560084863, 0.13002033209}

Out[436]= {{0.0777960697813, 0.556265936296, 0.807681686834, 
  0.420186553434}, {0.527982841158, 0.555876762387, 0.981090055861, 
  0.1560084863}, {0.907596678267, 0.364390115148, 0.722450567448, 
  0.13002033209}}

Out[437]= {{1, 0., 0., 0.415692234671}, {0, 1, 0., 4.46163163168}, {0, 0, 1, -2.59261340179}}

Check that this agrees with LinearSolve (I assume that was the purpose).

rred[[All, -1]] == LinearSolve[aA, b]

Out[438]= True

My apologies, I forgot about the last column as the title mentioned 'numerical errors'. Thanks for completing the answer.

Wow, it works. Why? What's the magic behind that?