How to implement a matrix-free method on sparse matrix? (Win11 x86-64 cuda)

[kevlns]

I want to implement an implicit viscosity solver which involves a 3nx3n matrix A, the linear system to calculated is (I - dtA)v = v0, where v and v0 are 3nx1 matrix. Actually A consists of 3x3 blocks, v and v0 consist of 3x1 blocks.