Add support for Apple's Accelerate sparse matrix solvers
What does this implement/fix?
This merge request implements support for the sparse matrix solvers found in Apple's Accelerate framework. Wrappers around the following solvers are provided:
- AccelerateLLT: Cholesky (LL^T) factorization
- AccelerateLDLT: Default LDL^T factorization (Currently an alias to AccelerateLDLTTPP)
- AccelerateLDLTUnpivoted: Cholesky-like LDL^T with only 1x1 pivots and no pivoting
- AccelerateLDLTSBK: LDL^T with Supernode Bunch-Kaufman and static pivoting
- AccelerateLDLTTPP: LDL^T with full threshold partial pivoting
- AccelerateQR: QR factorization
- AccelerateCholeskyAtA: QR factorization without storing Q (equivalent to A^TA = R^T R)
Performance
Notes:
- Default solver settings are used in each test.
- Eigen QR did not finish in a reasonable amount of time.
- Tests run on an Apple M1 Devkit (4e4p cores)
Solving system of size: 1946848x1946848
| Solver | Analysis (ms) | Factorization (ms) | Solve (ms) | Total (ms) |
|---|---|---|---|---|
| AccelerateLDLT | 722 | 9861 | 169 | 10753 |
| AccelerateLLT | 744 | 774 | 242 | 1760 |
| Eigen LDLT | 1376 | 45685 | 343 | 47405 |
| Eigen LLT | 1366 | 45716 | 346 | 47429 |
| Cholmod Simplicial LDLT | 12134 | 10171 | 79 | 22385 |
| Cholmod Supernodal LLT | 13138 | 864 | 110 | 14114 |
Solving system of size: 57504x57504:
| Solver | Analysis (ms) | Factorization (ms) | Solve (ms) | Total (ms) |
|---|---|---|---|---|
| AccelerateQR | 221 | 2939 | 217 | 3378 |
| SPQR | 69 | 961 | 4371 | |
| Eigen QR | 64 | > 60000 |
Additional information
In order for Accelerate to work, you need to provide which triangle you wish to use (for LDLT and LLT) and the type of matrix that you are factorizing. Eg symmetric, ordinary or triangular. The way that I'm determining this currently is by the exploiting the UpLo template argument.
The question I have is what would be a sensible UpLo default? If you try and run Accelerate's LDLT solver and m_sparseKind is not SparseSymmetric, it will assert and crash the program. In the way that I'm currently handling it, you would therefore need to do: AccelerateLDLT<SparseMatrix<double>, Symmetric | Upper> LDLT; This doesn't exactly follow what the other solvers do however. We could always assume that whatever you're passing into the solver is positive definite, which is what the PARDISO wrapper seems to do. This would work for LDLT and LLT, but not necessarily for QR.
QR is a bit tricky as your matrix may not be upper or lower triangular, and may also not be symmetric - ie it's "ordinary" as Apple calls it. However UpLoType does not have an appropriate entry for this case, so you must pass 0 in order to set m_sparseKind to SparseOrdinary. Eg: AccelerateQR<SparseMatrix<double>, 0> QR; You can see an example of this in the test program.