@Override
protected boolean convergedI(double r, Vector x)
        throws IterativeSolverNotConvergedException {
    // Store initial residual
    if (isFirst())
        initR = r;

    // Check for convergence
    if (r < Math.max(rtol * (normA * x.norm(normType) + normb), atol))
        return true;

    // Check for divergence
    if (r > dtol * initR)
        throw new IterativeSolverNotConvergedException(
                NotConvergedException.Reason.Divergence, this);
    if (iter >= maxIter)
        throw new IterativeSolverNotConvergedException(
                NotConvergedException.Reason.Iterations, this);
    if (Double.isNaN(r))
        throw new IterativeSolverNotConvergedException(
                NotConvergedException.Reason.Divergence, this);

    // Neither convergence nor divergence
    return false;
}
 

@Override
protected boolean convergedI(double r)
        throws IterativeSolverNotConvergedException {
    // Store initial residual
    if (isFirst())
        initR = r;

    // Check for convergence
    if (r < Math.max(rtol * initR, atol))
        return true;

    // Check for divergence
    if (r > dtol * initR)
        throw new IterativeSolverNotConvergedException(
                NotConvergedException.Reason.Divergence, this);
    if (iter >= maxIter)
        throw new IterativeSolverNotConvergedException(
                NotConvergedException.Reason.Iterations, this);
    if (Double.isNaN(r))
        throw new IterativeSolverNotConvergedException(
                NotConvergedException.Reason.Divergence, this);

    // Neither convergence nor divergence
    return false;
}
 

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

    double rho_1 = 1, rho_2 = 1, alpha = 1, beta = 1;

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

    for (iter.setFirst(); !iter.converged(r, x); iter.next()) {
        M.apply(r, z);
        M.transApply(rtilde, ztilde);
        rho_1 = z.dot(rtilde);

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

        if (iter.isFirst()) {
            p.set(z);
            ptilde.set(ztilde);
        } else {
            beta = rho_1 / rho_2;
            p.scale(beta).add(z);
            ptilde.scale(beta).add(ztilde);
        }

        A.mult(p, q);
        A.transMult(ptilde, qtilde);

        alpha = rho_1 / ptilde.dot(q);
        x.add(alpha, p);
        r.add(-alpha, q);
        rtilde.add(-alpha, qtilde);

        rho_2 = rho_1;
    }

    return x;
}
 

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

    double rho_1 = 1, rho_2 = 1, alpha = 1, beta = 1, omega = 1;

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

    for (iter.setFirst(); !iter.converged(r, x); iter.next()) {
        rho_1 = rtilde.dot(r);

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

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

        if (iter.isFirst())
            p.set(r);
        else {
            beta = (rho_1 / rho_2) * (alpha / omega);

            // temp = p - omega * v
            temp.set(-omega, v).add(p);

            // p = r + beta * temp = r + beta * (p - omega * v)
            p.set(r).add(beta, temp);
        }

        M.apply(p, phat);
        A.mult(phat, v);
        alpha = rho_1 / rtilde.dot(v);
        s.set(r).add(-alpha, v);

        x.add(alpha, phat);
        if (iter.converged(s, x))
            return x;

        M.apply(s, shat);
        A.mult(shat, t);
        omega = t.dot(s) / t.dot(t);
        x.add(omega, shat);
        r.set(s).add(-omega, t);

        rho_2 = rho_1;
    }

    return x;
}
 

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

    double rho_1 = 0, rho_2 = 0, alpha = 0, beta = 0;

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

    for (iter.setFirst(); !iter.converged(r, x); iter.next()) {
        rho_1 = rtilde.dot(r);

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

        if (iter.isFirst()) {
            u.set(r);
            p.set(u);
        } else {
            beta = rho_1 / rho_2;
            u.set(r).add(beta, q);
            sum.set(q).add(beta, p);
            p.set(u).add(beta, sum);
        }

        M.apply(p, phat);
        A.mult(phat, vhat);
        alpha = rho_1 / rtilde.dot(vhat);
        q.set(-alpha, vhat).add(u);

        M.apply(sum.set(u).add(q), uhat);
        x.add(alpha, uhat);
        A.mult(uhat, qhat);
        r.add(-alpha, qhat);

        rho_2 = rho_1;
    }

    return x;
}
 

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;
}