no.uib.cipr.matrix.NotConvergedException Java Examples

The following examples show how to use no.uib.cipr.matrix.NotConvergedException. You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may check out the related API usage on the sidebar.
Example #1
Source File: MatrixIterationMonitor.java    From matrix-toolkits-java with GNU Lesser General Public License v3.0 6 votes vote down vote up 
@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;
}
Example #2
Source File: DefaultIterationMonitor.java    From matrix-toolkits-java with GNU Lesser General Public License v3.0 6 votes vote down vote up 
@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;
}
Example #3
Source File: BiCG.java    From matrix-toolkits-java with GNU Lesser General Public License v3.0 5 votes vote down vote up 
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;
}
Example #4
Source File: BiCGstab.java    From matrix-toolkits-java with GNU Lesser General Public License v3.0 4 votes vote down vote up 
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;
}
Example #5
Source File: CGS.java    From matrix-toolkits-java with GNU Lesser General Public License v3.0 4 votes vote down vote up 
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;
}
Example #6
Source File: QMR.java    From matrix-toolkits-java with GNU Lesser General Public License v3.0 4 votes vote down vote up 
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;
}