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

Gordon,

Replace -Inverse[A].b with -b.Inverse[A] and then use Chop as suggested by Guillermo.

A = Table[RandomReal[], {i, 3}, {j, 3}] ;
b = Table[RandomReal[], {i, 3}] ;
M = Transpose[Join[A, {b}]] ;
ret = M.Join[
    Transpose[Join[Inverse[A], {-b.Inverse[A]}]], {{0, 0, 0, 1}}] ;
ret // MatrixForm
ret // Chop // MatrixForm

I.M.

It's so good to learn ArrayFlatten. Helps a lot.

Thank you