I am not sure exactly what you are trying to achieve, but some points:
If you want to deal with the matrices yourself you can do so in Python using petsc4py or even numpy, and I am confident that there is some C++ equivalent (using PETSc directly?). However...
I am not sure if I understand the 'to sum or merge them into one matrix' part. Instead of complex eigenvalue problem you should have two real ones, one for the real part and one for the imaginary part, where there might be some coupled terms present (i.e. dependency on the imaginary solution in the evp of the real part or vice versa). The best way to deal with this (Python API, you'll have to look for the C++ equivalent)
V = FunctionSpace(mesh, 'CG' 1)
Vcomplex = V * V
ureal, uimag = TrialFunction(Vcomplex)
vreal, vimag = TestFunction(Vcomplex)
Then you can write the weak form with both real and imaginary parts together, where you use vreal
for the real terms and vimag
for the imaginary ones. When you assemble, it will automatically result in 1 matrix.