Make sure the sparse solver is used when appropriate

[casella]

Please check the attached test package. It contains a battery of scaled tests involving sparse initial equations systems of increasing size. There are on average about three variables in each equation, a degree of sparsity that is commonly found in object-oriented models. Tearing is deactivated, to inquire the behaviour of the naked numerical solver.

These are the reported times spent on initialization with increasing size, with the default option and with simulation flags -nls=kinsol, -nlsLS=klu

N	T default	T flags
10	0.0004	0.00018
20	0.0015	0.00023
40	0.004	0.00028
80	0.047	0.00049
160	0.085	0.00076
320	0.52	0.0015
640	6.7	0.0028

The default values for the automatic activation of the sparse nonlinear solver are -nlssMinSize=10001 and -nlssMaxDensity=0.2. In all the cases reported in the second column, when the defaults are not changed, the used solver is the built-in dense Newton/homotopy solver.

It is very clearly visible as the solution times grows as O(N³) in the default case, because of the computational complexity of dense linear solvers. On the other hand, the time spent by the sparse solver grows as O(N).

For this very sparse system (the density is about 3/N for large enough N) the sparse solver systematically outperforms the dense one. Until about 100 equations both are very fast anyway, but the time spent by the dense solver goes quickly through the roof as the size goes about that, making it very inefficient to solve systems above 300 equations. This is way below the current default limit.

We should probably try out sparse systems with higher density and sparsity patterns other than tri-diagonal, but I have no time to do that right now. My current experience with ScalableTestGrid systems, where initial equations mostly involve two or three variables each, confirm the findings on the attached simple test case.

I would then drastically reduce the default value of -nlssMinSize, and also of -lssMinSize, since the performance penalty of the nonlinear solver is entirely due to the poor performance of the linear solver employed by the Newton algorithm, so there is no reason to treat linear and nonlinear systems in different ways.

Based on the results of this test, we could take a default value which is as low as 10. Recall that we also have a maximum limit of density -nlssMaxDensity=0.2, so systems of slightly higher size, say 20 or 30 equations, but also with density higher than 0.2, e.g. because of tearing, would still be handled by the dense solver, which would probably be the most appropriate in that case.

[casella]

See ​PR 7083.

[casella]

OK, things are a bit more complicated than I thought.

On one hand, there is probably an incentive at using the builtin/dense Newton/homotopy solver for initialization, as it has a number of heuristics and tricks to handle convergence failures that Kinsol does not have. On the other hand, in case there are large algebraic loops in the simulation equations, the best solver is IDA in daeMode (see also #6340), so this option won't be used during simulation, as daeMode can only work with IDA/Kinsol and needs to use a sparse solver for systems with more than a few dozens variables (i.e. practically all meaningful Modelica models).

So, it makes sense to increase the default minimium size for the nonlinear solver to a larger size; even though the dense solver is less efficient, it is only used once at initialiation, and there is a premium at having increased robustness. On the other hand, initialization should not take an unreasonable amount of time. Given the results shown in the table, I think a reasonable compromise is to set nlssMinSize=100 by default.

There is no such issue in the case of linear solvers; in fact, if we turn tearing off by default (see #6196 and ​PR 6969) it will be essential not to use the dense solver inappropriately. So, I'll leave the limit to 10 there. In fact, we could even set it to a lower value.

I updated ​PR 7083 accordingly.

[casella]

There is another issue, though, which is demonstrated in the updated test package, DefaultSolverReducedMinSize package. Even if I set -nlssMinSize=10, the builtin Newton/homotopy solver is still used; apparently, since it only works with a dense linear solver, it ignores that flag.

In order to get the sparse nonlinear solver to do its job, I need to set -nls="kinsol", and at that point, apparently, it uses the sparse linear solver by default anyway.

This is all a bit messy, I guess we should clean this up for 1.17.0 if possible. We can discuss this during the next devmeeting.

[casella]

Retargeted to 1.18.0 because of 1.17.0 timed release.

[casella]

I just updated the attached sample library, which now includes tests for both linear and nonlinear systems.

I confirm that with the default settings, sparse solvers are never selected, even in cases where it would be absolutely obvious that it is advantageous to do so.

[casella]

I made some further enquiry with the attached test suite. For example, if I run Linear.DefaultSolver.Test_N_80, it uses Lapack, i.e. the default.

If I then add -lssMinSize=1 I get this interesting output

Maximum system size for using linear sparse solver changed to 1

initialize linear system solvers
1 linear systems

Using sparse solver for linear system 0,
because density of 0.037 remains under threshold of 0.200 and size of 80 exceeds threshold of 1.

The maximum density and the minimal system size for using sparse solvers can be specified
using the runtime flags '<-lssMaxDensity=value>' and '<-lssMinSize=value>'.

Start solving Linear System 161 (size 80) at time 0 with Klu Solver

So it seems that the logic to select the sparse solver is:

if size > lssMinSize and density < lssMaxDensity then use sparse else use dense

I don't find this very useful. I mean, lssMinSize is set by default to 4001. I don't think it is reasonable to use a dense solver when you have 4000 unknowns: the jacobian would take 16 M elements, that is to say 128 MB of RAM. If the density is below the threshold, we should use a sparse solver regardless of the size. If the size is above the threshold, we should use a sparse solver anyway, to try minimizing the memory footprint.

So, to my understanding, the selection logic should be:

if size > lssMinSize or density < lssMaxDensity then use sparse else use dense

Ditto for nonlinear solvers

if size > nlssMinSize or density < nlssMaxDensity then use sparse else use dense

We could then trim the defaults a bit based on some experiments, but we should change the logic right away. The current one often caused some of my models to become impossibly slow, and besides, without giving any clue so as to the cause. They were simply trying to factorize a dense matrix with some million zeros in it.

[casella]

Done in ​PR 7695.

@phannebohm, can you please port this to v.18 maintenance as well.

[casella]

Ticket retargeted after milestone closed

[casella]

Issue eventually resolved in 1.19.0