finally I've managed to create some working code. http://ubuntuone.com/6nFP4CkylZXjFRaswmoP6S
is a small video which shows the results, a little assembly created with constraints. Some constraint types are still missing and it's not stable in all cases, but I think this is a big step.
As mentioned before I've rewritten the solver to get it working, changed it to the Levenberg Marquat algorithm. But this was still not enough, so i decided to change the solver to a two-step solution: first calculate all rotational constraints, after that all translational. That works pretty nice, but actually the old solver could do that as well. I've found my implementation to be more reliable (especially with the translational constraints), but maybe we need a deep comparison here.
@logari: you've mentioned stability issues with the single-subsystem solver before, maybe you could try my implementation on a few test cases with known issues? I've you're interested I could create a patch!
My next steps will be to make this more stable, clean up my code and again make it more stable