课程中有mass spring 牛顿法求解的代码， 所以我目前在尝试写FEM牛顿法求解的代码。
问题是我把整个hessian矩阵算出来了，但是Hessian 矩阵太大必须使用稀疏存储。
而ti.linalg.SparseMatrixBuilder 只能在cpu上运行，Sparse Spatial Data Structures可以在gpu上使用，但是并没有对应的矩阵计算模块。我测试了一下，同样的张量计算求Hessian矩阵的值，gpu的速度是cpu的大约30倍。请问如果想要用gpu加速稀疏Hessian矩阵的计算应该如何处理？
可以考虑用 matrix free 的方法，不用存下整个矩阵，而是对于每个四面体计算 hessian，在做矩阵乘法时累加所有四面体对于结果的贡献。可以参考 taichi/implicit_fem.py at master · taichi-dev/taichi · GitHub ，也即 ti example implicit_fem。
谢谢你的回答。matrix-free方法是个不错的方案，不过我之后考虑复现ipc，matrix-free的公式太难推广到其他情况了。目前我看到比较好的解决方法是直接不构建hessian矩阵。用文献Parallel Iterative Solvers for Real-time Elastic Deformations中的Preconditioning Descent Method by Jacobi+Chebyshev方法。
另外吐个槽，implicit_fem的代码cg迭代部分没有注释是真的太难看懂了。