I think came upon a bug in the mkl_?omatadd routine, especially the complex versions. If both matrices A and B are conjugate transposed, the matrix A is only conjugated. I configured the zomatadd.c example file that intel provided in order to show that this is the case:
As you can see the results are the same, which is incorrect. The second computation gives the wrong answer.
There are actually two more cases I found, where the result is wrong:
The example file is also attached for you to test it.