Support sparse matrix and sparse linear solvers

[FantasyVR]

Concisely describe the proposed feature Many problems in computer graphics (especially physical simulation) need to solve sparse linear systems. While some of these problems can be addressed using matrix-free solvers (e.g., conjugate gradients), many need explicitly constructing a sparse matrix, in standard formats such as CSR and COO.

Currently, Taichi does not have first-class support for these matrices. Note that the sparse computation system in Taichi is designed for spatial sparsity instead of matrix sparsity.

Describe the solution you'd like (if any)

• [x] Use SparseMatrixBuilder to create triplets of sparse matrices.
• [x] SparseMatrix is built from SparseMatrixBuilder
  • The users could configure the storage format (CSR, COO, etc.).
  • The users could configure the data type of sparse matrix element.
• [x] Support basic sparse matrix operations such as +, -, *, @, transpose, etc.
• [x] Support sparse direct solvers such as LLT, LDLT, LU, etc.
  • The users could configure the data type, ordering of solvers.
• [x] Performance benchmark (Scipy, cuSolver)
• [x] cuSPARSE, cuSolver function loader
• [x] create a sparse matrix using cuSparse
• [ ] Support GPU-based basic sparse matrix operations such as +, -, *, @, transpose, etc.
• [x] COO format sort before transforming to CSR format, reference: here.
• [ ] Datatype configuration of GPU sparse matrix/sparse solver.
• [ ] Performance benchmark of GPU/CPU sparse matrix/solver.

---

Update:

• [ ] Support more common used packages

  • SuiteSparse, Pardiso

• [ ] Remove SparseMatrixBuilder and build SparseMatrix using taichi ndarray

[mzy2240]

Will you consider supporting sparse matrix solver like Suitesparse?

[zishengye]

Will you consider an inclusion of some engineering level packages like PETSc and Trilinos. These packages not only support different types of matrix formats but also a lot of preconditioners.

[FantasyVR]

Will you consider supporting sparse matrix solvers like Suitesparse?

@mzy2240 Thanks for your suggestion. We just wrap Eigen for now. SuiteSparse is in our plan. If you have interest,
Welcome to join us and support more features.

[FantasyVR]

Will you consider the inclusion of some engineering level packages like PETSc and Trilinos. These packages not only support different types of matrix formats but also a lot of preconditioners.

Our plan is to design a general sparse linear system solver. But now, we intend to implement a basic and useable version based on Eigen. More packages like SuitSparse, Pardiso, CuSolver et al. are in our plan. And we welcome you to join us during the process.

[zishengye]

Maybe KokkosKernel could help your work.

[mzy2240]

Do you consider supporting conversion from and to scipy.sparse matrix? Scipy is probably the most used library when dealing with sparse matrix in python.

[Hanke98]

Hi, there! I am now working on implementing the basic operation operations such as +, -, *, @, etc. on GPU. However, I found that the + and - operations are not common cases according to this post.

Is it OK to use cuSparse Routine csrgeam, which performs C=α∗A+β∗B, to implement the + operation?

[FantasyVR]

Hi @Hanke98. Exactly. I thinkcsrgeam could be used to implement +, -, scalar *. I just implemented a not finished one. It could be as a reference.