这是一些用于Java的BLAS库的很好的摘要:performance-of-java-matrix-math-
libraries。您还可以在Java-Matrix-Benchmark上查看其中许多库的基准。
但是,根据我的经验,这些库中的大多数似乎都不适合解决大型稀疏矩阵。在我的情况下,我所做的是通过JNI
使用Eigen实现解决方案。
Eigen对其线性求解器进行了很好的讨论,其中包括一个在CHOLMOD上的讨论。
以我为例,使用Eigen的求解器通过JNI生成的8860x8860稀疏矩阵比并行柯尔特快20倍,比我自己的密集求解器快10倍。更重要的是,它看起来像按比例缩放,
n^2而不是按比例缩放,
n^3并且它使用的内存比密集型求解器要少得多(我用尽了按比例扩大内存)。
实际上,Eigen有一个Java封装器,称为JEigen,它使用JNI。但是,它没有实现稀疏矩阵求解,因此不会包装所有内容。
我最初使用JNA,但对开销不满意。Wikipedia
很好地说明了如何使用JNI。一旦编写了函数声明并与之一起编译,就
javac可以使用
javah它们创建C
++的头文件。
例如
//Cholesky.javapackage cfd.optimisation;//ri, ci, v : matrix row indices, column indices, and values//y = Ax where A is a nxn matrix with nnz non-zero valuespublic 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 DoubleMatrix2Dint 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);


