哪位大佬能介绍一下matrix-free方法？

综合讨论区
#1

搞流体模拟的隐式粘度求解器，直接求解线性系统是真算不动啊。网上没搜到很直观的关于matrix-free方法的介绍啊555（具体如 matrix-free conjugate gradient。有好一些的参考资料吗？）

#2

Matrix-free 的求解器可以参考一下 Scipy 的文档：scipy.sparse.linalg.LinearOperator — SciPy v1.10.1 Manual

目前 Taichi 里也有了一个 matrix-free 的 CG 求解器：[lang] Implement experimental CG(Conjugate Gradient) solver in Taichi-lang by houkensjtu · Pull Request #7690 · taichi-dev/taichi · GitHub

使用方法暂时可以参考测试脚本：taichi/test_taichi_cg.py at master · taichi-dev/taichi · GitHub

简单说就是你可以定义一个 LinearOperator 来代替矩阵 A ，然后传入 CG 求解器来求解方程组。有问题的话欢迎继续讨论！

1 Like
#3

非常感谢，我先学习学习~

#4

这个测试脚本已经失效，请问还有别的参考嘛

#5

谢谢指出！我来整理一下最新的版本稍后发在 repo 里～

#6

最近的测试脚本在这里：https://github.com/taichi-dev/taichi/blob/master/tests/python/test_matrixfree_cg.py

下面是一个直接可运行的版本：

import taichi as ti
from taichi.linalg import taichi_cg_solver, LinearOperator # Use taichi_cg_solver in 1.6.0
# from taichi.linalg import MatrixFreeCG, LinearOperator # Use MatrixFreeCG from 1.7.0 (not released yet)

ti.init(arch=ti.cpu)

n = 4
x = ti.field(dtype=ti.f32, shape=(n, n))
Ax = ti.field(dtype=ti.f32, shape=(n, n))
b = ti.field(dtype=ti.f32, shape=(n, n))

@ti.kernel
def init():
    for i in ti.grouped(b):
        b[i] = 1.0
        x[i] = 0.0

@ti.kernel
def compute(v:ti.template(), mv:ti.template()):
    for i in ti.grouped(v):
        mv[i] = 2 * v[i]

init()        
A = LinearOperator(compute)
print(b)
taichi_cg_solver(A, b, x, maxiter=10*n*n, tol=1e-18, quiet=True)
# MatrixFreeCG(A, b, x, maxiter=10*n*n, tol=1e-18, quiet=True)  # Use MatrixFreeCG if you're on master branch
print(x)

注意如果你是 <=1.6.0 的版本，那么需要使用 taichi_cg_solver 这个调用名。如果是源码编译的，调用名是 MatrixFreeCG。另外，1.6.0 版本还有一个 bug 导致只能求解 2 维问题，新的 master 上已经修复了不过还没反映在 pip 版本里。