/** */ private static void axpy(Double a, double[] x, double[] y) { f2jBlas.daxpy(x.length, a, x, 1, y, 1); }
/** * @param M tall, skinny matrix * @return MT * M as a dense lower-triangular matrix, represented in packed row-major form. */ public static double[] transposeTimesSelf(Collection<float[]> M) { if (M == null || M.isEmpty()) { return null; } int features = M.iterator().next().length; double[] result = new double[features * (features + 1) / 2]; double[] vectorD = new double[features]; for (float[] vector : M) { // Unfortunately need to copy into a double[] for (int i = 0; i < vector.length; i++) { vectorD[i] = vector[i]; } blas.dspr("L", features, 1.0, vectorD, 1, result); } return result; }
/** {@inheritDoc} */ @Override protected double[] iter(double bnorm, double[] target) { double[] res = dataset.computeWithCtx((ctx, data) -> { if (data.getFeatures() == null) return null; int cols = data.getFeatures().length / data.getRows(); BLAS.getInstance().dscal(ctx.getU().length, 1 / bnorm, ctx.getU(), 1); double[] v = new double[cols]; BLAS.getInstance().dgemv("T", data.getRows(), cols, 1.0, data.getFeatures(), Math.max(1, data.getRows()), ctx.getU(), 1, 0, v, 1); return v; }, (a, b) -> { if (a == null) return b; else if (b == null) return a; else { BLAS.getInstance().daxpy(a.length, 1.0, a, 1, b, 1); return b; } }); BLAS.getInstance().daxpy(res.length, 1.0, res, 1, target, 1); return target; }
/** {@inheritDoc} */ @Override protected double beta(double[] x, double alfa, double beta) { return dataset.computeWithCtx((ctx, data) -> { if (data.getFeatures() == null) return null; int cols = data.getFeatures().length / data.getRows(); BLAS.getInstance().dgemv("N", data.getRows(), cols, alfa, data.getFeatures(), Math.max(1, data.getRows()), x, 1, beta, ctx.getU(), 1); return BLAS.getInstance().dnrm2(ctx.getU().length, ctx.getU(), 1); }, (a, b) -> a == null ? b : b == null ? a : Math.sqrt(a * a + b * b)); }
alfa = blas.dnrm2(v.length, v, 1); blas.dscal(v.length, 1 / alfa, v, 1); anorm = Math.sqrt(Math.pow(anorm, 2) + Math.pow(alfa, 2) + Math.pow(beta, 2) + Math.pow(damp, 2)); blas.dscal(v.length, -beta, v, 1); alfa = blas.dnrm2(v.length, v, 1); blas.dscal(v.length, 1 / alfa, v, 1); double t2 = -theta / rho; blas.dcopy(w.length, w, 1, dk, 1); blas.dscal(dk.length, 1 / rho, dk, 1); blas.daxpy(w.length, t1, w, 1, x, 1); blas.dscal(w.length, t2, w, 1); blas.daxpy(w.length, 1, v, 1, w, 1); ddnorm = ddnorm + Math.pow(blas.dnrm2(dk.length, dk, 1), 2); blas.daxpy(var.length, 1.0, pow2(dk), 1, var, 1);
/** * Returns true if BLAS is available in any form. This should always return * true because we've included JBLAS. * * @return true if BLAS is available in any form. */ public static boolean blasAvailable() { return (BLAS.getInstance() != null); }
/** {@inheritDoc} */ @Override protected double bnorm() { return dataset.computeWithCtx((ctx, data) -> { ctx.setU(Arrays.copyOf(data.getLabels(), data.getLabels().length)); return BLAS.getInstance().dnrm2(data.getLabels().length, data.getLabels(), 1); }, (a, b) -> a == null ? b : b == null ? a : Math.sqrt(a * a + b * b)); }
/** * Sets up an eigenvalue decomposition for symmetrical, tridiagonal * matrices. Uses a low default tolerance criteria * * @param n * Size of the matrix * @param vectors * True to compute the eigenvectors, false for just the * eigenvalues */ public SymmTridiagEVD(int n, boolean vectors) { this(n, vectors, LAPACK.getInstance().dlamch("Safe minimum")); }
/** * For the moment we have no flags indicating if matrix is transposed or not. Therefore all dgemm parameters for * transposition are equal to 'N'. */ public static void gemm(double alpha, Matrix a, Matrix b, double beta, Matrix c) { if (alpha == 0.0 && beta == 1.0) return; else if (alpha == 0.0) scal(c, beta); else { double[] fA = a.getStorage().data(); double[] fB = b.getStorage().data(); double[] fC = c.getStorage().data(); assert fA != null; nativeBlas.dgemm("N", "N", a.rowSize(), b.columnSize(), a.columnSize(), alpha, fA, a.rowSize(), fB, b.rowSize(), beta, fC, c.rowSize()); if (c instanceof SparseMatrix) MatrixUtil.unflatten(fC, c); } }
/** * Performs in-place multiplication of vector x by a real scalar a. (x = a * x) * * @param a Scalar a. * @param x Vector x. **/ public static void scal(Double a, Vector x) { if (x.isArrayBased()) f2jBlas.dscal(x.size(), a, x.getStorage().data(), 1); else if (x instanceof SparseVector) { Set<Integer> indexes = ((SparseVector)x).indexes(); for (Integer i : indexes) x.compute(i, (ind, v) -> v * a); } else throw new IllegalArgumentException(); }
/** * y := alpha * A * x + beta * y. * * @param alpha Alpha. * @param a Matrix a. * @param x Vector x. * @param beta Beta. * @param y Vector y. */ public static void gemv(double alpha, Matrix a, Vector x, double beta, Vector y) { checkCardinality(a, x); if (a.rowSize() != y.size()) throw new CardinalityException(a.columnSize(), y.size()); if (alpha == 0.0 && beta == 1.0) return; if (alpha == 0.0) { scal(y, beta); return; } double[] fA = a.getStorage().data(); double[] fX = x.getStorage().data(); double[] fY = y.getStorage().data(); nativeBlas.dgemv("N", a.rowSize(), a.columnSize(), alpha, fA, a.rowSize(), fX, 1, beta, fY, 1); if (y instanceof SparseVector) y.assign(fY); }
/** TODO: IGNTIE-5770, add description for a */ static void syr(Double alpha, DenseVector x, DenseMatrix a) { int nA = a.rowSize(); int mA = a.columnSize(); nativeBlas.dsyr("U", x.size(), alpha, x.getStorage().data(), 1, a.getStorage().data(), nA); // Fill lower triangular part of A int i = 0; while (i < mA) { int j = i + 1; while (j < nA) { a.setX(j, i, a.getX(i, j)); j++; } i++; } }
@Override public void dtptri(java.lang.String uplo, java.lang.String diag, int n, double[] ap, int _ap_offset, org.netlib.util.intW info) { dtptri_offsets(uplo, diag, n, ap, _ap_offset, info); }
@Override public void sggrqf(int m, int p, int n, float[] a, int _a_offset, int lda, float[] taua, int _taua_offset, float[] b, int _b_offset, int ldb, float[] taub, int _taub_offset, float[] work, int _work_offset, int lwork, org.netlib.util.intW info) { sggrqf_offsets(m, p, n, a, _a_offset, lda, taua, _taua_offset, b, _b_offset, ldb, taub, _taub_offset, work, _work_offset, lwork, info); }
/** * Returns true if BLAS is available in any form. This should always return * true because we've included JBLAS. * * @return true if BLAS is available in any form. */ public static boolean blasAvailable() { return (BLAS.getInstance() != null); }
/** * Sets up an eigenvalue decomposition for symmetrical, tridiagonal * matrices. Uses a low default tolerance criteria * * @param n * Size of the matrix * @param vectors * True to compute the eigenvectors, false for just the * eigenvalues */ public SymmTridiagEVD(int n, boolean vectors) { this(n, vectors, LAPACK.getInstance().dlamch("Safe minimum")); }
/** * Adds alpha * v * v.t to a matrix in-place. This is the same as BLAS's ?SPR. * * @param u the upper triangular part of the matrix in a [[DenseVector]](column major) */ public static void spr(Double alpha, DenseVector v, DenseVector u) { nativeBlas.dspr("U", v.size(), alpha, v.getStorage().data(), 1, u.getStorage().data()); }
/** * Calculates standard deviation by all columns. * * @return Standard deviations. */ public double[] std() { double[] mean = mean(); ValueWithCount<double[]> res = delegate.compute(data -> { double[] features = data.getFeatures(); int rows = data.getRows(); int cols = features.length / rows; double[] y = new double[cols]; for (int col = 0; col < cols; col++) for (int j = col * rows; j < (col + 1) * rows; j++) y[col] += Math.pow(features[j] - mean[col], 2); return new ValueWithCount<>(y, rows); }, (a, b) -> a == null ? b : b == null ? a : new ValueWithCount<>(sum(a.val, b.val), a.cnt + b.cnt)); if (res != null) { blas.dscal(res.val.length, 1.0 / res.cnt, res.val, 1); for (int i = 0; i < res.val.length; i++) res.val[i] = Math.sqrt(res.val[i]); return res.val; } return null; }
/** * Returns true if BLAS is available in any form. This should always return * true because we've included JBLAS. * * @return true if BLAS is available in any form. */ public static boolean blasAvailable() { return (BLAS.getInstance() != null); }
/** * Calculates mean value by all columns. * * @return Mean values. */ public double[] mean() { ValueWithCount<double[]> res = delegate.compute((data, partIdx) -> { double[] features = data.getFeatures(); int rows = data.getRows(); int cols = features.length / rows; double[] y = new double[cols]; for (int col = 0; col < cols; col++) for (int j = col * rows; j < (col + 1) * rows; j++) y[col] += features[j]; return new ValueWithCount<>(y, rows); }, (a, b) -> a == null ? b : b == null ? a : new ValueWithCount<>(sum(a.val, b.val), a.cnt + b.cnt)); if (res != null) { blas.dscal(res.val.length, 1.0 / res.cnt, res.val, 1); return res.val; } return null; }