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

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

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

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

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

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

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

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

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

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