Switch default solver to Richardson iteration

[pierrelux]

Our current default solver amounts to Richardson iteration with scaling factor "omega" equals to 1. Divergence is likely to occur if the spectral radius of df/dx in f(x,theta) = x is too close to 1. Note that under the assumption that we have a fixed point process and that "f" is a contraction mapping, there shouldn't be an issue with the existence of the inverse jacobian in the IFT. However, things get more complicated when we depart from the pure FP perspective and try to support both FP and root finding under the same API.