I want to solve a decently sized DAE system (up to 1e4 equations) with IDA which is too slow (and of course memory consuming) using dense linear algebra operations on the Jacobian. Therefore I either need to use the direct sparse solver module or the iterative methods of SUNDIALS. I successfully used the former with KINSOL in a C++ project with 1e5 (and more) equations but I have to rely on Python for my current project and the direct solvers are currently not usable in Assimulo (as far as I know?). Unfortunenately I can't get the GMRES solver to work either, probably because of the lack of preconditioning which is currently also not available for IDA.
At the moment I'm using IDA because of a singular mass matrix. However the system can be approximated using a modified non-singular mass matrix and therefore my next step will be trying to solve the modified system with CVODE. Do you have any other ideas regarding solving ODE/DAE systems with sparse Jacobian matrices in Python?