Java Code Examples for no.uib.cipr.matrix.Vector#add()

The following examples show how to use no.uib.cipr.matrix.Vector#add() . 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: LinkedSparseMatrix.java    From matrix-toolkits-java with GNU Lesser General Public License v3.0 5 votes vote down vote up 
@Override
public Vector multAdd(double alpha, Vector x, Vector y) {
    checkMultAdd(x, y);
    if (alpha == 0)
        return y;
    Node node = links.head;
    while (node != null) {
        y.add(node.row, alpha * node.val * x.get(node.col));
        node = node.rowTail;
    }
    return y;
}
Example 2
Source File: LinkedSparseMatrix.java    From matrix-toolkits-java with GNU Lesser General Public License v3.0 5 votes vote down vote up 
@Override
public Vector transMultAdd(double alpha, Vector x, Vector y) {
    checkTransMultAdd(x, y);
    if (alpha == 0)
        return y;
    Node node = links.head;
    while (node != null) {
        y.add(node.col, alpha * node.val * x.get(node.row));
        node = node.colTail;
    }
    return y;
}
Example 3
Source File: FlexCompRowMatrix.java    From matrix-toolkits-java with GNU Lesser General Public License v3.0 5 votes vote down vote up 
@Override
public Vector multAdd(double alpha, Vector x, Vector y) {
    checkMultAdd(x, y);

    for (int i = 0; i < numRows; ++i)
        y.add(i, alpha * rowD[i].dot(x));

    return y;
}
Example 4
Source File: IR.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);

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

    for (iter.setFirst(); !iter.converged(r, x); iter.next()) {
        M.apply(r, z);
        x.add(z);
        A.multAdd(-1, x, r.set(b));
    }

    return x;
}
Example 5
Source File: CG.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 alpha = 0, beta = 0, rho = 0, rho_1 = 0;

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

    for (iter.setFirst(); !iter.converged(r, x); iter.next()) {
        M.apply(r, z);
        rho = r.dot(z);

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

        A.mult(p, q);
        alpha = rho / p.dot(q);

        x.add(alpha, p);
        r.add(-alpha, q);

        rho_1 = rho;
    }

    return x;
}
Example 6
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 7
Source File: Chebyshev.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 alpha = 0, beta = 0, c = 0, d = 0;

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

    c = (eigmax - eigmin) / 2.0;
    d = (eigmax + eigmin) / 2.0;

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

        if (iter.isFirst()) {
            p.set(z);
            alpha = 2.0 / d;
        } else {
            beta = (alpha * c) / 2.0;
            beta *= beta;
            alpha = 1.0 / (d - beta);
            p.scale(beta).add(z);
        }

        A.mult(p, q);
        x.add(alpha, p);
        r.add(-alpha, q);
    }

    return x;
}
Example 8
Source File: FlexCompColMatrix.java    From matrix-toolkits-java with GNU Lesser General Public License v3.0 5 votes vote down vote up 
@Override
public Vector transMultAdd(final double alpha, final Vector x,
        final Vector y) {
    checkTransMultAdd(x, y);

    for (int i = 0; i < numColumns; ++i)
        y.add(i, alpha * colD[i].dot(x));

    return y;
}
Example 9
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 10
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 11
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;
}
Example 12
Source File: GMRES.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);

    A.multAdd(-1, x, u.set(b));
    M.apply(u, r);
    double normr = r.norm(Norm.Two);
    M.apply(b, u);

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

        v[0].set(1 / normr, r);
        s.zero().set(0, normr);
        int i = 0;

        // Inner iteration
        for (; i < restart && !iter.converged(Math.abs(s.get(i))); i++, iter
                .next()) {
            A.mult(v[i], u);
            M.apply(u, w);

            for (int k = 0; k <= i; k++) {
                H.set(k, i, w.dot(v[k]));
                w.add(-H.get(k, i), v[k]);
            }
            H.set(i + 1, i, w.norm(Norm.Two));
            v[i + 1].set(1. / H.get(i + 1, i), w);

            // QR factorization of H using Givens rotations
            for (int k = 0; k < i; ++k)
                rotation[k].apply(H, i, k, k + 1);

            rotation[i] = new GivensRotation(H.get(i, i), H.get(i + 1, i));
            rotation[i].apply(H, i, i, i + 1);
            rotation[i].apply(s, i, i + 1);
        }

        // Update solution in current subspace
        new UpperTriangDenseMatrix(H, i, false).solve(s, s);
        for (int j = 0; j < i; j++)
            x.add(s.get(j), v[j]);

        A.multAdd(-1, x, u.set(b));
        M.apply(u, r);
        normr = r.norm(Norm.Two);
    }

    return x;
}