FeynCalc Commutator of Matrices

If I understand the question correctly, what you need is Inner[Dot, x, y]

Here is an example of working with Dirac matrices by representing them as 2x2 matrices containing Pauli matrices. From here on you should be able to cook up the code you need for your problem.

gamma[0] = {{1, 0}, {0, -1}}
gamma[i_] := {{0, CSI[i]}, {-CSI[i], 0}}

blockMatrixProduct[x_] := x;
blockMatrixProduct[x_, y_] := Inner[Dot, x, y];
blockMatrixProduct[x_, y_, z__] := 
blockMatrixProduct[x, blockMatrixProduct[y, z]];

blockMatrixProduct[gamma[i], gamma[j], gamma[i]]
DotSimplify[%]

$$ left( begin{array}{cc} 0.left(overline{sigma }^j.left(-overline{sigma }^iright)+0.0right)+overline{sigma }^i.left(0.left(-overline{sigma }^iright)+left(-overline{sigma }^jright).0right) & 0.left(0.overline{sigma }^i+overline{sigma }^j.0right)+overline{sigma }^i.left(left(-overline{sigma }^jright).overline{sigma }^i+0.0right) 0.left(0.left(-overline{sigma }^iright)+left(-overline{sigma }^jright).0right)+left(-overline{sigma }^iright).left(overline{sigma }^j.left(-overline{sigma }^iright)+0.0right) & 0.left(left(-overline{sigma }^jright).overline{sigma }^i+0.0right)+left(-overline{sigma }^iright).left(0.overline{sigma }^i+overline{sigma }^j.0right) end{array} right) $$

$$ left( begin{array}{cc} 0 & -overline{sigma }^i.overline{sigma }^j.overline{sigma }^i overline{sigma }^i.overline{sigma }^j.overline{sigma }^i & 0 end{array} right) $$