[REQ] Direct Sparse Cholesky Factorization and Implicit Euler

[yuhan-zh]

Description

Hi! I am inquiring if Warp has any plans to add a direct sparse solver for large Cholesky Factorization, potentially by leveraging the cuDSS library. A direct benefit of this would be the ability to support Implicit Euler solvers.

Context

I am aware that Warp is not primarily designed for simulations using fully implicit solvers with large global systems. And there are already XPBD and VBD which are implicit method that can converge to the solution of implicit Euler with newton's method theoretically(while they normally cannot). However, implicit solvers play a crucial role in physical simulation, and it would be very valuable to have this capability in Warp as an option or for testing/benchmarking purposes(for example, to compare with the result from XPBD or VBD).

Thank you!

[gdaviet]

Hi @yuhan-zh , user-friendly bindings to cuDSS is indeed something we want to have from warp.sparse eventually.

In the meantime, depending on your problem the iterative linear solvers from warp.optim.linear can be an alternative. There are also a few examples that demonstrate using them as part of implicit euler solvers, for instance https://github.com/NVIDIA/warp/blob/bf8830d11365ac5061c779c55b5278d733076d4f/warp/examples/fem/example_mixed_elasticity.py#L184

[yuhan-zh]

Hi @gdaviet, Thank you for your fast reply and suggestion! May I ask is there a rough timeline for adding the cuDSS bindings? I am mainly wondering if it's something that is going to be available in maybe a few months or probably more than a year.

Thank you!

[shi-eric]

I'll also mention that the nvmath-python library already exposes cuDSS bindings: https://docs.nvidia.com/cuda/nvmath-python/latest/bindings/cudss.html but I haven't tried them myself.

[gdaviet]

Thanks @shi-eric. I've asked cursor to generate an example using this API and it seems to work fine. We could potentially add a csr_to_cupy helper to streamline this, though the verbosity is really not too bad:

    # Create CuPy arrays using zero-copy interface
    csr_indptr = cp.asarray(matrix_bsr.offsets)
    csr_indices = cp.asarray(matrix_bsr.columns)
    # Reshape values from (nnz, 1, 1) to (nnz,) for standard CSR format
    csr_data = cp.asarray(matrix_bsr.values).reshape(-1)

    # Create CuPy CSR matrix
    A_cupy = cp.sparse.csr_matrix((csr_data, csr_indices, csr_indptr), shape=(N, N), dtype=cp.float32)

@momo-van for viz

Full example: example_warp_sparse_cudss.py

[yuhan-zh]

Thank you @gdaviet @shi-eric. This is a very useful information and a nice workaround. And thank you for the example.

Still, I hope there will be natively supported sparse solver with adjoint computation in Warp someday. Thanks again!