ЧОЛМОД на Java

Вот хороший обзор некоторых библиотек BLAS для Java: математические библиотеки. Вы также можете ознакомиться с тестами многих из этих библиотек на Java-Matrix-Benchmark.

Однако, по моему опыту, большинство этих библиотек не приспособлены для решения больших разреженных матриц. В моем случае я реализовал решение, используя Eigen через JNI.

У Eigen есть хорошее обсуждение своих линейных решателей, в том числе и CHOLMOD.

В моем случае для разреженной матрицы 8860x8860 использование решателя Эйгена через JNI было в 20 раз быстрее, чем параллельный кольт, и в 10 раз быстрее, чем мой собственный плотный решатель. Более важно то, что он масштабируется как n^2, а не n^3, и использует гораздо меньше памяти, чем мой плотный решатель (мне не хватило памяти при масштабировании).

На самом деле существует оболочка для Eigen с Java под названием JEigen, которая использует JNI. Однако в нем не реализовано решение разреженных матриц, поэтому он не охватывает все.

Сначала я использовал JNA, но не был доволен накладными расходами. В Википедии есть хороший пример использования JNI. После того, как вы написали объявления функций и скомпилировали их с помощью javac, вы используете javah для создания файла заголовка для C++.

Например для

//Cholesky.java
package cfd.optimisation;
//ri, ci, v : matrix row indices, column indices, and values
//y = Ax where A is a nxn matrix with nnz non-zero values
public class Cholesky {
    private static native void solve_eigenLDLTx(int[] ri, int[] ci, double[] v, double[] x, double[] y, int n, int nnz);
}

С помощью javah создается заголовочный файл cfd_optimization_Cholesky.h с объявлением

JNIEXPORT void JNICALL Java_cfd_optimisation_Cholesky_solve_1eigenLDLTx
        (JNIEnv *, jclass, jintArray, jintArray, jdoubleArray, jdoubleArray, jdoubleArray, jint, jint); 

И вот как я реализовал решатель

JNIEXPORT void JNICALL Java_cfd_optimisation_Cholesky_solve_1eigenLDLTx(JNIEnv *env, jclass obj, jintArray arrri, jintArray arrci, jdoubleArray arrv, jdoubleArray arrx, jdoubleArray arry, jint jn, jint jnnz) {
    int n = jn;
    int *ri = (int*)env->GetPrimitiveArrayCritical(arrri, 0);
    int *ci = (int*)env->GetPrimitiveArrayCritical(arrci, 0);
    double *v = (double*)env->GetPrimitiveArrayCritical(arrv, 0);
    int nnz = jnnz;

    double *x = (double*)env->GetPrimitiveArrayCritical(arrx, 0);
    double *y = (double*)env->GetPrimitiveArrayCritical(arry, 0);

    Eigen::SparseMatrix<double> A = colt2eigen(ri, ci, v, nnz, n);
    //Eigen::MappedSparseMatrix<double> A(n, n, nnz, ri, ci, v);

    Eigen::VectorXd a(n), b(n);
    for (int i = 0; i < n; i++) a(i) = x[i];
    //a = Eigen::Map<Eigen::VectorXd>(x, n).cast<double>(); 
    Eigen::SimplicialCholesky<Eigen::SparseMatrix<double> > solver;
    solver.setMode(Eigen::SimplicialCholeskyLDLT);
    b = solver.compute(A).solve(a);
    for (int i = 0; i < n; i++) y[i] = b(i);
    env->ReleasePrimitiveArrayCritical(arrri, ri, 0);
    env->ReleasePrimitiveArrayCritical(arrci, ci, 0);
    env->ReleasePrimitiveArrayCritical(arrv, v, 0);
    env->ReleasePrimitiveArrayCritical(arrx, x, 0);
    env->ReleasePrimitiveArrayCritical(arry, y, 0);
}

Функция colt2eigen создает разреженную матрицу из двух целочисленных массивов, содержащих индексы строк и столбцов, и двойной массив значений.

Eigen::SparseMatrix<double> colt2eigen(int *ri, int *ci, double* v, int nnz, int n) {
    std::vector<Eigen::Triplet<double>> tripletList;
    for (int i = 0; i < nnz; i++) { 
        tripletList.push_back(Eigen::Triplet<double>(ri[i], ci[i], v[i]));  
    }
    Eigen::SparseMatrix<double> m(n, n);
    m.setFromTriplets(tripletList.begin(), tripletList.end());
    return m;
}

Одной из сложных частей было получение этих массивов от Java и Colt. Для этого я сделал это

//y = A x: x and y are double[] arrays and A is DoubleMatrix2D
int nnz = A.cardinality();
DoubleArrayList v = new DoubleArrayList(nnz);
IntArrayList ci = new IntArrayList(nnz);
IntArrayList ri = new IntArrayList(nnz);

A.forEachNonZero((row, column, value) -> {
    v.add(value); ci.add(column); ri.add(row); return value;}
);

Cholesky.solve_eigenLDLTx(ri.elements(), ci.elements(), v.elements(), x, y, n, nnz);

Я не работал ни с одним из этих инструментов, но я подозреваю, что вас поражает тот факт, что Java не использует встроенную инструкцию процессора с плавающей точкой квадратного корня в некоторых версиях/на некоторых платформах.

См.: Где я могу найти источник код для функции квадратного корня Java?

Если абсолютная точность не требуется, вы можете попробовать переключить одну из приведенных выше реализаций на использование приближения квадратного корня (см. Быстрый sqrt в Java за счет точности для предложений), что должно быть несколько быстрее.