/** {@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; }
/** * 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); } }
/** * @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 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)); }
/** {@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)); }
@Override public Matrix multAdd(double alpha, Matrix B, Matrix C) { if (!(B instanceof DenseMatrix) || !(C instanceof DenseMatrix)) return super.multAdd(alpha, B, C); checkMultAdd(B, C); double[] Bd = ((DenseMatrix) B).getData(), Cd = ((DenseMatrix) C) .getData(); BLAS.getInstance().dgemm(Transpose.NoTranspose.netlib(), Transpose.NoTranspose.netlib(), C.numRows(), C.numColumns(), numColumns, alpha, data, Math.max(1, numRows), Bd, Math.max(1, B.numRows()), 1, Cd, Math.max(1, C.numRows())); return C; }
@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(); BLAS.getInstance().dgemv(Transpose.NoTranspose.netlib(), numRows, numColumns, alpha, data, Math.max(numRows, 1), xd, 1, 1, yd, 1); return y; }
/** * Calculate the 2 norm of the ndarray * * @param A * @return */ public static double nrm2(IComplexNDArray A) { if (A.data().dataType() == DataBuffer.Type.FLOAT) { float s = BLAS.getInstance().snrm2(A.length(), A.data().asFloat(), 2, 1); return s; } else if (A.data().dataType() == DataBuffer.Type.DOUBLE) { double s = BLAS.getInstance().dnrm2(A.length(), A.data().asDouble(), 2, 1); return s; } throw new IllegalArgumentException("Illegal data type"); }
/** * Copy x to y * * @param x the origin * @param y the destination */ public static void copy(IComplexNDArray x, IComplexNDArray y) { DataTypeValidation.assertSameDataType(x, y); if (x.data().dataType() == DataBuffer.Type.FLOAT) BLAS.getInstance().scopy( x.length(), x.data().asFloat(), x.majorStride(), y.data().asFloat(), y.majorStride()); else BLAS.getInstance().dcopy( x.length(), x.data().asDouble(), x.majorStride(), y.data().asDouble(), y.majorStride()); }
@Override public Matrix rank1(double alpha, Vector x, Vector y) { if (!(x instanceof DenseVector) || !(y instanceof DenseVector)) return super.rank1(alpha, x, y); checkRank1(x, y); double[] xd = ((DenseVector) x).getData(), yd = ((DenseVector) y) .getData(); BLAS.getInstance().dger(numRows, numColumns, alpha, xd, 1, yd, 1, data, Math.max(1, numRows)); return this; }
@Override public void daxpy(int n,double alpha, DataBuffer x, int offsetX, int incrX, DataBuffer y, int offsetY, int incrY ){ BLAS.getInstance().daxpy(n, alpha, getDoubleData(x), offsetX, incrX, getDoubleData(y), offsetY, incrY); }
@Override public Matrix rank1(double alpha, Vector x, Vector y) { if (x != y) throw new IllegalArgumentException("x != y"); if (!(x instanceof DenseVector)) return super.rank1(alpha, x, y); checkRank1(x, y); double[] xd = ((DenseVector) x).getData(); BLAS.getInstance().dspr(uplo.netlib(), numRows, alpha, xd, 1, data); return this; }
@Override protected void dscal(int N, double alpha, INDArray X, int incX) { double[] data = getDoubleData(X); BLAS.getInstance().dscal(N, alpha, data, BlasBufferUtil.getBlasOffset(X), incX); setData(data,X); }
@Override protected double ddot( int N, DataBuffer X, int offsetX, int incX, DataBuffer Y, int offsetY, int incY){ return BLAS.getInstance().ddot(N, getDoubleData(X), offsetX, incX, getDoubleData(Y), offsetY, incY); }
DataTypeValidation.assertSameDataType(x, y); if (x.data().dataType() == DataBuffer.Type.FLOAT) { double ret = BLAS.getInstance().sdot( x.length(), x.data().asFloat(), double ret = BLAS.getInstance().ddot( x.length(), x.data().asDouble(),
@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(); BLAS.getInstance().dgbmv(Transpose.Transpose.netlib(), numRows, numColumns, kl, ku, alpha, data, kl + ku + 1, xd, 1, 1, yd, 1); 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(), yd = ((DenseVector) y) .getData(); BLAS.getInstance().dsbmv(uplo.netlib(), numRows, kd, alpha, data, kd + 1, xd, 1, 1, yd, 1); return y; }
@Override protected void dcopy(int n, DataBuffer x, int offsetX, int incrX, DataBuffer y, int offsetY, int incrY ){ BLAS.getInstance().dcopy(n, getDoubleData(x), offsetX, incrX, getDoubleData(y), offsetY, incrY); }
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); }