no.uib.cipr.matrix.Vector Java Examples

@Override
public Vector multAdd(double alpha, Vector x, Vector y) {
    if (!(x instanceof DenseVector) || !(y instanceof DenseVector))
        return super.multAdd(alpha, x, y);

    checkMultAdd(x, y);

    double[] xd = ((DenseVector) x).getData(), yd = ((DenseVector) y)
            .getData();

    // y = 1/alpha * y
    y.scale(1 / alpha);

    // y = A*x + y
    for (int i = 0; i < numColumns; ++i)
        for (int j = columnPointer[i]; j < columnPointer[i + 1]; ++j)
            yd[rowIndex[j]] += data[j] * xd[i];

    // y = alpha*y = alpha*A*x + y
    return y.scale(alpha);
}
 

@Override
public Vector transMult(Vector x, Vector y) {
    if (!(x instanceof DenseVector) || !(y instanceof DenseVector))
        return super.transMult(x, y);

    checkTransMultAdd(x, y);

    double[] xd = ((DenseVector) x).getData();
    double[] yd = ((DenseVector) y).getData();

    y.zero();

    for (int i = 0; i < numRows; ++i)
        for (int j = rowPointer[i]; j < rowPointer[i + 1]; ++j)
            yd[columnIndex[j]] += data[j] * xd[i];

    return y;
}
 

/**
 * Constructor for SparseVector, and copies the contents from the supplied
 * vector.
 * 
 * @param x
 *            Vector to copy from
 * @param deep
 *            True if a deep copy is to be made. If the copy is shallow,
 *            <code>x</code> must be a <code>SparseVector</code>
 */
public SparseVector(Vector x, boolean deep) {
    super(x);

    if (deep) {
        int nz = Matrices.cardinality(x);
        data = new double[nz];
        index = new int[nz];
        set(x);
    } else {
        SparseVector xs = (SparseVector) x;
        data = xs.getData();
        index = xs.getIndex();
        used = xs.getUsed();
    }
}
 

@Override
public Vector mult(Vector x, Vector y) {
    if (!(x instanceof DenseVector) || !(y instanceof DenseVector))
        return super.mult(x, y);

    checkMultAdd(x, y);

    double[] xd = ((DenseVector) x).getData();
    double[] yd = ((DenseVector) y).getData();

    y.zero();

    for (int i = 0; i < ind.length; ++i) {
        int row = ind[i] < 0 ? -ind[i] : 0;
        int column = ind[i] > 0 ? ind[i] : 0;
        double[] locDiag = diag[i];
        for (int j = 0; j < locDiag.length; ++j, ++row, ++column)
            yd[row] += locDiag[j] * xd[column];
    }

    return y;
}
 

@Override
public Vector multAdd(double alpha, Vector x, Vector y) {
    if (!(x instanceof DenseVector) || !(y instanceof DenseVector))
        return super.multAdd(alpha, x, y);

    checkMultAdd(x, y);

    double[] xd = ((DenseVector) x).getData();
    double[] yd = ((DenseVector) y).getData();

    for (int i = 0; i < ind.length; ++i) {
        int row = ind[i] < 0 ? -ind[i] : 0;
        int column = ind[i] > 0 ? ind[i] : 0;
        double[] locDiag = diag[i];
        for (int j = 0; j < locDiag.length; ++j, ++row, ++column)
            yd[row] += alpha * locDiag[j] * xd[column];
    }

    return y;
}
 

@Override
public Vector transMultAdd(double alpha, Vector x, Vector y) {
    if (!(x instanceof DenseVector) || !(y instanceof DenseVector))
        return super.transMultAdd(alpha, x, y);

    checkTransMultAdd(x, y);

    double[] xd = ((DenseVector) x).getData();
    double[] yd = ((DenseVector) y).getData();

    for (int i = 0; i < ind.length; ++i) {
        int row = ind[i] < 0 ? -ind[i] : 0;
        int column = ind[i] > 0 ? ind[i] : 0;
        double[] locDiag = diag[i];
        for (int j = 0; j < locDiag.length; ++j, ++row, ++column)
            yd[column] += alpha * locDiag[j] * xd[row];
    }

    return y;
}
 

@Override
public Vector multAdd(double alpha, Vector x, Vector y) {
    if (!(x instanceof DenseVector) || !(y instanceof DenseVector))
        return super.multAdd(alpha, x, y);

    checkMultAdd(x, y);

    double[] xd = ((DenseVector) x).getData();
    double[] yd = ((DenseVector) y).getData();

    // y = 1/alpha * y
    y.scale(1. / alpha);

    // y = A*x + y
    for (int i = 0; i < numColumns; ++i) {
        SparseVector v = colD[i];
        int[] index = v.getIndex();
        double[] data = v.getData();
        int length = v.getUsed();
        for (int j = 0; j < length; ++j)
            yd[index[j]] += data[j] * xd[i];
    }

    // y = alpha*y = alpha * A'x + y
    return y.scale(alpha);
}
 

/**
 * weights(cat1) ... weights(cat2) ... ... bias(cat1) bias(cat2)
 */
public double[] getWeightVector() {
	final double[] result = new double[(totalFeatures + 1) * lexicalCategories.size()];
	int index = 0;
	for (final Vector vector : weightMatrixRows) {
		for (int i = 0; i < vector.size(); i++) {
			result[index] = vector.get(i);
			index++;
		}
	}

	for (int i = 0; i < bias.size(); i++) {
		result[index] = bias.get(i);
		index++;
	}

	return result;
}
 

@Override
public Vector transMultAdd(double alpha, Vector x, Vector y) {
    if (!(x instanceof DenseVector) || !(y instanceof DenseVector))
        return super.transMultAdd(alpha, x, y);

    checkTransMultAdd(x, y);

    double[] xd = ((DenseVector) x).getData();
    double[] yd = ((DenseVector) y).getData();

    // y = 1/alpha * y
    y.scale(1. / alpha);

    // y = A'x + y
    for (int i = 0; i < numRows; ++i)
        for (int j = rowPointer[i]; j < rowPointer[i + 1]; ++j)
            yd[columnIndex[j]] += data[j] * xd[i];

    // y = alpha*y = alpha*A'x + y
    return y.scale(alpha);
}
 

@Override
public Vector multAdd(double alpha, Vector x, Vector y) {
    if (!(x instanceof DenseVector) || !(y instanceof DenseVector))
        return super.multAdd(alpha, x, y);

    checkMultAdd(x, y);

    double[] xd = ((DenseVector) x).getData();
    double[] yd = ((DenseVector) y).getData();

    for (int i = 0; i < numRows; ++i) {
        double dot = 0;
        for (int j = rowPointer[i]; j < rowPointer[i + 1]; ++j)
            dot += data[j] * xd[columnIndex[j]];
        yd[i] += alpha * dot;
    }

    return y;
}
 

/**
 * Constructor for QMR. Uses the given vector as template for creating
 * scratch vectors. Typically, the solution or the right hand side vector
 * can be passed, and the template is not modified. Allows setting different
 * right and left preconditioners
 * 
 * @param template
 *            Vector to use as template for the work vectors needed in the
 *            solution process
 * @param M1
 *            Left preconditioner
 * @param M2
 *            Right preconditioner
 */
public QMR(Vector template, Preconditioner M1, Preconditioner M2) {
    this.M1 = M1;
    this.M2 = M2;
    r = template.copy();
    y = template.copy();
    z = template.copy();
    v = template.copy();
    w = template.copy();
    p = template.copy();
    q = template.copy();
    d = template.copy();
    s = template.copy();
    v_tld = template.copy();
    w_tld = template.copy();
    y_tld = template.copy();
    z_tld = template.copy();
    p_tld = template.copy();
}
 

@Override
public Vector transSolve(Vector b, Vector x) {
    if (!(x instanceof DenseVector))
        return super.transSolve(b, x);

    x.set(b);

    double[] xd = ((DenseVector) x).getData();

    for (int i = numRows - 1; i >= 0; --i) {

        // Get row i
        SparseVector row = LU.getRow(i);
        int[] index = row.getIndex();
        double[] data = row.getData();

        // At this stage, x[i] is known, so move it over to the right
        // hand side for the remaining equations
        for (int j = 0; j < diagind[i]; ++j)
            xd[index[j]] -= data[j] * xd[i];

    }

    return x;
}
 

@Override
public Vector transSolve(Vector b, Vector x) {
    if (!(x instanceof DenseVector))
        return super.transSolve(b, x);

    x.set(b);

    double[] xd = ((DenseVector) x).getData();

    for (int i = 0; i < numRows; ++i) {

        // Solve for the current entry
        xd[i] /= data[diagind[i]];

        // Move this known solution over to the right hand side for the
        // remaining equations
        for (int j = diagind[i] + 1; j < rowptr[i + 1]; ++j)
            xd[colind[j]] -= data[j] * xd[i];
    }

    return x;
}
 

@Override
public Vector solve(Vector b, Vector x) {
    if (!(b instanceof DenseVector) || !(x instanceof DenseVector))
        return super.solve(b, x);

    double[] bd = ((DenseVector) b).getData();
    double[] xd = ((DenseVector) x).getData();

    for (int i = numRows - 1; i >= 0; --i) {

        // xi = (bi - sum[j>i] Uij * xj) / Uii
        double sum = 0;
        for (int j = diagind[i] + 1; j < rowptr[i + 1]; ++j)
            sum += data[j] * xd[colind[j]];

        xd[i] = (bd[i] - sum) / data[diagind[i]];
    }

    return x;
}
 

@Override
public Vector solve(Vector b, Vector x) {
    if (!(b instanceof DenseVector) || !(x instanceof DenseVector))
        return super.solve(b, x);

    double[] bd = ((DenseVector) b).getData();
    double[] xd = ((DenseVector) x).getData();

    for (int i = 0; i < numRows; ++i) {

        // xi = bi - sum[j<i] Lij * xj
        double sum = 0;
        for (int j = rowptr[i]; j < diagind[i]; ++j)
            sum += data[j] * xd[colind[j]];

        xd[i] = bd[i] - sum;
    }

    return x;
}
 

@Override
public Vector transSolve(Vector b, Vector x) {
    if (!(x instanceof DenseVector))
        return super.transSolve(b, x);

    x.set(b);

    double[] xd = ((DenseVector) x).getData();

    for (int i = numRows - 1; i >= 0; --i)

        // At this stage, x[i] is known, so move it over to the right hand
        // side for the remaining equations
        for (int j = rowptr[i]; j < diagind[i]; ++j)
            xd[colind[j]] -= data[j] * xd[i];

    return x;
}
 

@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
public Vector transMult(Vector x, Vector y) {
    if (!(x instanceof DenseVector) || !(y instanceof DenseVector))
        return super.transMult(x, y);

    checkTransMultAdd(x, y);

    double[] xd = ((DenseVector) x).getData();
    double[] yd = ((DenseVector) y).getData();

    for (int i = 0; i < numColumns; ++i) {
        double dot = 0;
        for (int j = columnPointer[i]; j < columnPointer[i + 1]; ++j)
            dot += data[j] * xd[rowIndex[j]];
        yd[i] = dot;
    }

    return y;
}
 

@Override
public Vector transMultAdd(double alpha, Vector x, Vector y) {
    if (!(x instanceof DenseVector) || !(y instanceof DenseVector))
        return super.transMultAdd(alpha, x, y);

    checkTransMultAdd(x, y);

    double[] xd = ((DenseVector) x).getData(), yd = ((DenseVector) y)
            .getData();

    // y = 1/alpha * y
    y.scale(1. / alpha);

    // y = A'x + y
    for (int i = 0; i < numRows; ++i) {
        SparseVector v = rowD[i];
        int[] index = v.getIndex();
        double[] data = v.getData();
        int length = v.getUsed();
        for (int j = 0; j < length; ++j)
            yd[index[j]] += data[j] * xd[i];
    }

    // y = alpha*y = alpha * A'x + y
    return y.scale(alpha);
}
 

@Override
public Vector solve(Vector b, Vector x) {
    if (!(b instanceof DenseVector) || !(x instanceof DenseVector))
        return super.solve(b, x);

    double[] bd = ((DenseVector) b).getData();
    double[] xd = ((DenseVector) x).getData();

    for (int i = 0; i < numRows; ++i) {

        // Get row i
        SparseVector row = LU.getRow(i);
        int[] index = row.getIndex();
        double[] data = row.getData();

        // xi = bi - sum[j<i] Lij * xj
        double sum = 0;
        for (int j = 0; j < diagind[i]; ++j)
            sum += data[j] * xd[index[j]];

        xd[i] = bd[i] - sum;
    }

    return x;
}
 

@Override
public Vector solve(Vector b, Vector x) {
    if (!(b instanceof DenseVector) || !(x instanceof DenseVector))
        return super.solve(b, x);

    double[] bd = ((DenseVector) b).getData();
    double[] xd = ((DenseVector) x).getData();

    for (int i = numRows - 1; i >= 0; --i) {

        // Get row i
        SparseVector row = LU.getRow(i);
        int[] index = row.getIndex();
        int used = row.getUsed();
        double[] data = row.getData();

        // xi = (bi - sum[j>i] Uij * xj) / Uii
        double sum = 0;
        for (int j = diagind[i] + 1; j < used; ++j)
            sum += data[j] * xd[index[j]];

        xd[i] = (bd[i] - sum) / data[diagind[i]];
    }

    return x;
}
 

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

@Override
void testFactorization(Matrix A, Vector x) {
    Vector b = A.mult(x, x.copy());

    ILU ilu = new ILU(new CompRowMatrix(A));
    ilu.setMatrix(A);
    ilu.apply(b, x);

    Vector r = A.multAdd(-1, x, b.copy());

    assertEquals(0, r.norm(Vector.Norm.TwoRobust), 1e-5);
}
 

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

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 transApply(Vector b, Vector x) {
    // U'y = b, y = U'\b
    U.transSolve(b, y);

    // L'x = U'\b = y
    return L.transSolve(y, x);
}
 

public Vector apply(Vector b, Vector x) {
    // Ly = b, y = L\b
    L.solve(b, y);

    // Ux = L\b = y
    return U.solve(y, 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 apply(Vector b, Vector x) {
    if (!(x instanceof DenseVector) || !(b instanceof DenseVector))
        throw new IllegalArgumentException("Vector must be DenseVectors");

    double[] xd = ((DenseVector) x).getData();
    double[] bd = ((DenseVector) b).getData();

    for (int i = 0; i < invdiag.length; ++i)
        xd[i] = bd[i] * invdiag[i];

    return x;
}
 

/**
 * Constructor for CGS. Uses the given vector as template for creating
 * scratch vectors. Typically, the solution or the right hand side vector
 * can be passed, and the template is not modified
 * 
 * @param template
 *            Vector to use as template for the work vectors needed in the
 *            solution process
 */
public CGS(Vector template) {
    p = template.copy();
    q = template.copy();
    u = template.copy();
    phat = template.copy();
    qhat = template.copy();
    vhat = template.copy();
    uhat = template.copy();
    sum = template.copy();
    r = template.copy();
    rtilde = template.copy();
}