public Vector solve(Matrix A, Vector b, Vector x)
            throws IterativeSolverNotConvergedException {
        checkSizes(A, b, x);

        double rho = 0, rho_1 = 0, xi = 0, gamma = 1., gamma_1 = 0, theta = 0, theta_1 = 0, eta = -1., delta = 0, ep = 0, beta = 0;

        A.multAdd(-1, x, r.set(b));

        v_tld.set(r);
        M1.apply(v_tld, y);
        rho = y.norm(Norm.Two);

        w_tld.set(r);
        M2.transApply(w_tld, z);
        xi = z.norm(Norm.Two);

        for (iter.setFirst(); !iter.converged(r, x); iter.next()) {

            if (rho == 0)
                throw new IterativeSolverNotConvergedException(
                        NotConvergedException.Reason.Breakdown, "rho", iter);

            if (xi == 0)
                throw new IterativeSolverNotConvergedException(
                        NotConvergedException.Reason.Breakdown, "xi", iter);

            v.set(1 / rho, v_tld);
            y.scale(1 / rho);
            w.set(1 / xi, w_tld);
            z.scale(1 / xi);

            delta = z.dot(y);

            if (delta == 0)
                throw new IterativeSolverNotConvergedException(
                        NotConvergedException.Reason.Breakdown, "delta", iter);

            M2.apply(y, y_tld);
            M1.transApply(z, z_tld);

            if (iter.isFirst()) {
                p.set(y_tld);
                q.set(z_tld);
            } else {
                p.scale(-xi * delta / ep).add(y_tld);
                q.scale(-rho * delta / ep).add(z_tld);
            }

            A.mult(p, p_tld);

            ep = q.dot(p_tld);

            if (ep == 0)
                throw new IterativeSolverNotConvergedException(
                        NotConvergedException.Reason.Breakdown, "ep", iter);

            beta = ep / delta;

            if (beta == 0)
                throw new IterativeSolverNotConvergedException(
                        NotConvergedException.Reason.Breakdown, "beta", iter);

            v_tld.set(-beta, v).add(p_tld);
            M1.apply(v_tld, y);
            rho_1 = rho;
            rho = y.norm(Norm.Two);

            A.transMultAdd(q, w_tld.set(-beta, w));
            M2.transApply(w_tld, z);
            xi = z.norm(Norm.Two);

            gamma_1 = gamma;
            theta_1 = theta;
            theta = rho / (gamma_1 * beta);
            gamma = 1 / Math.sqrt(1 + theta * theta);

            if (gamma == 0)
                throw new IterativeSolverNotConvergedException(
                        NotConvergedException.Reason.Breakdown, "gamma", iter);

            eta = -eta * rho_1 * gamma * gamma / (beta * gamma_1 * gamma_1);

            if (iter.isFirst()) {
                d.set(eta, p);
                s.set(eta, p_tld);
            } else {
                double val = theta_1 * theta_1 * gamma * gamma;
                d.scale(val).add(eta, p);
                s.scale(val).add(eta, p_tld);
            }

            x.add(d);
            r.add(-1, s);
        }

        return x;
    }