tools for sparse matrices

There are a few special operations that use sparse matrices. In particular, for solving tri-diagonal equations

Some of the operations in the MatrixOP suite, may be of use to you.

Hello Ken,

I looked into this issue recently and (unless I am missing something) I got the impression that the code that handles sparse matrices in e.g., Intel MKL, does not provide a very compelling argument for implementation.  I am therefore curious about your application and wonder if you could help me understand which aspect of sparse matrix handling is important for you. 

As usual, I'd appreciate receiving feedback from any customer that uses or needs to use sparse matrix operations.

A.G.

support@wavemetrics.com

Apologies for the slow response - one of the key people working on this problem was out of town...

Our application for handling sparse matrices is in solving atmospheric chemistry models in which there are up to 15000 reactions which form a stiff system of equations. One such model is called MCM.

http://mcm.leeds.ac.uk/MCM/

This can be run efficiently in MATLAB apparently, which uses a library similar to CVODE. From what I understand, there are more options in the CVODE library which are specifically designed to deal with sparse matrices. In MATLAB, a model run might take 5 minutes or so, but in Igor, with more than about 1000 reactions, it was taking many hours to run. One key to MATLAB seems to be the 'sparse' command. (I don't know MATLAB myself, but one of our colleagues uses it for this purpose.)

There are several research groups involved in atmospheric chemistry that we work with who would find such tools to be very useful.

Thanks for your response.

I'm curious what format is used to store the sparse matrices.  Is it the trivial full representation or are the non-zero elements sufficiently sparse that it makes more sense to represent them as a triplet {row,col,value}.

AG

I think MATLAB uses a triplet.

Here's an example of using sparse().

https://www.mathworks.com/help/matlab/math/constructing-sparse-matrices…

Ken

I'd like to re-start discussion on this issue. This factor of x50 slower performance for Igor vs. Matlab for chemical kinetics solvers for very large systems (e.g. 20,000 reactions) is now a big practical issue for the atmospheric chemistry community, which as you know has many Igor users. We have recently released a free Igor-based easy-to-use kinetic solver (KinSim v4) for teaching and research, which a lot of people are using now:

https://pubs.acs.org/doi/10.1021/acs.jchemed.9b00033

http://tinyurl.com/kinsim-release 

At least 15 papers have been published using KinSim so far. We have verified that KinSim runs as fast as the FACSIMILE commercial solver which some people in our community use. But the current performance advantage of Matlab is critical for many research applications. This will interest 100+ Igor users (over time, many more), including many of the Aerodyne instrument users. If the performance bottleneck of Igor is not resolved, then people will opt instead for using the free Matlab alternative, and once they start using Matlab, they may just stay with that for everything and abandon Igor. So I would urge the Igor folks to work on this sooner rather than later. 

Thanks,

-Jose

Hello Jose,

Thanks for the feedback.

Some sparse matrix support will be in IP9.  It might be useful to know specifically what operations are required in your application.  For example, alpha*MatA*vecX+beta*vecY for sparse matrix A and dense vectors X and Y and constants alpha and beta are already implemented.

A.G

It is an ODE solver for stiff systems, i.e. a version of your IntegrateODE that can operate on sparse stiff systems. I.e.

dXi/dt = Sum_j (kij * Xi * Xj)

where only a few terms are active. E.g. for a system with 20,000 species, the equations for most species may only have a few non-zero terms on the equation above.

Does that help? (if not we can dig further).

Sparse matrix operations are similar to standard matrix operations but they are specialized for sparse systems.  These could be used in other calculations but they are completely independent of the ODE solver.  I will let the person in charge of IntegrateODE respond to your query.

AG

Jose-

I take it that you are using the stiff system solver, that is IntegrateODE/M=3 to select the BDF method). But inside the derivative function you need to do some matrix computations on sparse matrices. So you are using the fastest solution method, and the speed is limited by the speed of using dense matrix computation?

- Yes, correct, we use /M=3

(I just talked to Matt Coggon of NOAA to clarify the details below in my mind, he may chip in further later)

- So our understanding is that the very slow part is inside IntegrateODE. I think there you use a Newton solver that needs to estimate the Jacobian by calling the derivative function multiple times. Our impression is that Igor has implemented an older version of the CVODE library (e.g. you cite the 1994 version of that code's User guide). Newer versions have implemented a sparse calculation. In particular, looking at the 2019 manual: 

https://computing.llnl.gov/sites/default/files/public/cv_guide.pdf (from https://computing.llnl.gov/projects/sundials/sundials-software)

on p38, there is a mention of how supplying and using the sparsity pattern of the matrix makes the solution much more efficient. Our understanding is that Matlab also has just implemented CVODE, but that they have a newer version that does allow taking advantage of sparsity and thus it runs much much faster for very large systems.

So to first order, our understanding is that Igor's IntegrateODE could be greatly sped up by implementing the sparse routine(s) within the updated CVODE. 

Perhaps in addition to that, sparse calculations within our derivative routine would help, but our understanding is that this is much less important. 

Thanks,

-Jose

Thanks for clarifying that. I was hoping I could make this AG's responsibility :) He does matrix computations (and is possibly implementing sparse matrix support for Igor 9).

I am, indeed, using a pretty old version of CVODE. I will have to look into updating it.

I'm not sure I have time to get it into Igor 9, since it seems like something that might require a fair bit of work.

Ok, please keep us posted as you learn more about this, and the potential timeline. 

I have filed an enhancement request ticket so I won't forget. I'm signed up for the Sundials mailing list, so maybe that will also help me remember.

Thanks. In terms of a timeline, can you estimate whether this requires a month or two, or 6 months, or a year etc? The earlier the better obviously as a lot of users would benefit, and they may go elsewhere otherwise.

It's very hard to say before I've seen the newer Sundials package. I don't even know if I will understand it :)

If I do this, I presume I can count on you for some testing?

Certainly

Jose-

My guess is that what you need is what's in CVODES. Could you look over the overview page and comment? It's here: https://computing.llnl.gov/projects/sundials/cvodes

The sensitivity and adjoint capabilities described there would be very nice. But what's urgent is just the faster stiff ODE solver that takes advantage of sparsity. Not sure how extensive the CVODES package is, but I am hoping that you could first implement the sparsity-using solver soon, and then the rest of the package could come at a later time.

Hi John - I agree with Jose regarding the ability to implement the sparsity pattern in CVODE or CVODES. It seems that capabilities in chapter 4 of the CVODES user guide would be most useful for the problem we are trying to attack. The implementation of a user-defined Jacobian (which utilizes sparsity patterns, page 77) would likely help the most.

Matt

I have an update on this problem.

I have downloaded the SUNDIALS 5 code, and have successfully built the libraries on my Macintosh. The changes to Igor's IntegrateODE code are substantial, but I have managed to get it to work on the ultra-simple first-order example from the manual: DisplayHelpTopic "A First-Order Equation"

There are error tolerance issues I still need to work out, but I don't think that will ultimately be too difficult.

I also need to work out the organization of our source code to make it easy to build for both Windows and Macintosh.

Once those steps are completed, I can start working on adding some of the newer, more sophisticated features of SUNDIALS.

I could use some help here, if you can, in understanding Chapter 2 of the SUNDIALS manual. In particular, it discusses various solvers. In particular, it talks about three "direct" solvers: dense, band and sparse. The word "sparse" caught my attention!

In the paragraph after the listing of solvers, it says: "For large stiff systems, where direct methods are often not feasible, the combination of a BDF in- tegrator and a preconditioned Krylov method yields..." This seems like it implies that the sparse solver cannot be used with the BDF integrator. Have I misunderstood something?

Hi John - My understanding is that the manual is describing systems with "really big" systems of equations, such as a 3D chemical transport modeling. These iterative methods are most often used in finite element analysis (such as was is used in fluid mechanic software packages, such as COMSOL). For the systems that we will be working with, the direct solvers should be sufficient since we are thinking about a single set of differential equations (this is what is used by the Matlab analog to the box modeling tool we are developing).

Matt

Hi John - Also, my interpretation of the manual is that they sparse solvers can be used for all systems of equations, it's just that for really big systems (such as for finite element analysis), one would need to use an iterative method combined with a direct solver (e.g. GMRES + BDF) rather than a direct solver by itself (i.e., BDF).

Thanks, Matt. I'm working on it. Right now I'm trying to figure out how to shoehorn the original error tolerance system into the new version of SUNDIALS.

Having successfully integrated the latest SUNDIALS release into Igor 9, I went to work on integrating sparse solvers. I presume that what you want would either SUNLinSol_KLU or SUNLinSol_SuperLUMT. It seems like KLU is more mature or perhaps user-friendly, but if that's true then SuperLUMT must be pretty terrible.

In any case, I downloaded the KLU package, or rather, the SuiteSparse package. There's a ton of code with lots of dependencies. I'm afraid that my conclusion is that I can't support this for Igor 9. I may come to the same conclusion in the future, too.

My apologies. Thanks for all your cooperation in looking into the matter.

FYI for people following the topic: I just talked to John on the phone.

- He will send us a beta of Igor 9 which has the new SunDials routines into IntegrateODE, so that we can test its speed against the old SunDials code in Igor 8

- More importantly, we are trying to work together to see if we can find a way to get the KLU code in a form that would make this possible.