Is it possible to recover *the* decomposition that's ultimately employed by mldivide (backslash)?
2 views (last 30 days)
If I issue:
Then I can see the decision made by mldivide on how to solve against A. In my case it chooses MA57:
sp\: bandwidth = 26633+1+26633.
sp\: is A diagonal? no.
sp\: is band density (0.00299541) > bandden (0.5) to try banded solver? no.
sp\: is A triangular? no.
sp\: is A morally triangular? no.
sp\: is A a candidate for Cholesky (symmetric, real positive or negative diagonal)? no.
sp\: But A is real symmetric; try MA57.
I would like to cache the decomposition (to conduct more solves against the same matrix later). If I call:
D = decomposition(A)
decomposition with properties:
MatrixSize: [27045 27045]
Sifting through the documentation I can match up ldl to MA57, so I thought this meant I'm using the same decomposition as above.
However a few things hint otherwise:
- \ is much faster
tic;D = decomposition(A);toc
tic;D = decomposition(A,'ldl');toc
on my MacPro with 1.5TB memory and 28 cores produces:
Elapsed time is 1.226560 seconds.
Elapsed time is 2.343294 seconds.
Elapsed time is 2.440572 seconds.
- solving with the decomposition gives me a warning that \ doesn't:
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.766761e-23.
> In decomposition/warnIfIllConditioned (line 714)
In \ (line 404)
- Do decomposition and \ make the same decisions on the choice of decomposition to use?
- If not, is there way to get at whatever decomposition \ is actually using? (in my case this "faster" MA57 decomposition)
I'm imagining/wish for something like:
[x,D] = A\b;
where D would contain the reusable decomposition created and used during \.
Here's A and b data
Christine Tobler on 25 Jun 2020
The same decomposition is used, the difference is in the estimation of the condition number. MA57 provides an estimate of the condition number which is dependent on the right-hand side vector used. Since this wasn't practical to use with the decomposition object, we call another estimator which unfortunately is more expensive. It's value will also not usually match the estimate from MA57 (which changes significantly based on the right-hand side vector), which would be why you're not getting that warning consistently.
You could turn off computing the condition number using decomposition(A, 'CheckCondition', false), this should show a noteable speedup in constructing the object.