以下是一些 Java BLAS 库的总结:java-matrix-math-libraries 的性能 https://stackoverflow.com/questions/529457/performance-of-java-matrix-math-libraries。您还可以在以下位置查看其中许多库的基准测试:Java 矩阵基准测试 https://lessthanoptimal.github.io/Java-Matrix-Benchmark/.
然而,根据我的经验,这些库中的大多数似乎并未针对解决大型稀疏矩阵进行调整。就我而言,我所做的是使用以下方法实现解决方案Eigen https://en.wikipedia.org/wiki/Eigen通过JNI。
Eigen对其线性求解器进行了很好的讨论 http://eigen.tuxfamily.org/dox/group__TopicSparseSystems.html其中包括乔莫德 (Cholmod) 上的一个。
对于我的 8860x8860 稀疏矩阵的情况,通过 JNI 使用 Eigen 求解器比并行 colt 快 20 倍,比我自己的密集求解器快 10 倍。更重要的是,它似乎可以扩展为n^2
而不是n^3
而且它使用的内存比我的密集解算器少得多(我在扩展时内存不足)。
实际上,Java 中有一个 Eigen 的包装器,称为JEigen https://github.com/hughperkins/jeigen它使用 JNI。然而,它没有实现稀疏矩阵求解,因此它没有包装所有内容。
我最初使用 JNA,但对开销并不满意。维基百科有关于如何使用 JNI 的一个很好的例子 https://en.wikipedia.org/wiki/Java_Native_Interface。一旦你编写了函数声明并编译它们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);
}
Using 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);