SparseCSCInterface sparspaklu! check sparsity pattern #27
Conversation
sparsepaklu! now has a 'reuse_symbolic' option:
sparspaklu!(lu::SparseSolver, m::SparseMatricCSC; reuse_symbolic=true)
If 'reuse_symbolic=true' (the default) then the sparsity pattern of 'm' is checked
and an error generated if it does not match that of the matrix used to create 'lu'
This matches the interface for the Julia SparseArrays stdlib UMFPACK lu!
Also some bug fixes:
- sparspaklu! always mutate and update 'lu'
- ErrorExceptions were created, but not actually thrown ! Fix by
replacing ErrorException() with error()
Change ErrorException(...) to error(...) so that an exception is thrown in SpkSparseSolver.solve! SpkSparseSpdSolver.solve!
sparspaklu!(lu::SparseSolver, m::SparseMatrixCSC; allow_pattern_change=true) If `allow_pattern_change = true` (the default) the sparse matrix `m` may have a nonzero pattern different to that of the matrix used to create the LU factorization `lu`, in which case the ordering and symbolic factorization are updated. If `allow_pattern_change = false` an error is thrown if the nonzero pattern of `m` is different to that of the matrix used to create the LU factorization `lu`. If `lu` was created with option `factorize = false` then `lu` is always updated from `m` and `allow_pattern_change` is ignored.
| @@ -1,5 +1,5 @@ | |||
| module SparseCSCInterface | |||
| using SparseArrays, LinearAlgebra | |||
| import SparseArrays, LinearAlgebra | |||
There was a problem hiding this comment.
Is an import necessary? Are we adding methods?
There was a problem hiding this comment.
I looked at this the other way around:
import is safer than an unqualified using, with the only difference being that it avoids importing many unknown and unused methods into the namespace with possibility of name clashes etc?
In this case, the only effect of this change was to enforce SparseArrays.nnz which looks like what was intended, and caught a couple of places where an unqualified nnz was inadvertently used which then was potentially confusing (well confused me anyway!).
Many thanks for fantastic work making this package available btw!
There was a problem hiding this comment.
In this case, the only effect of this change was to enforce
SparseArrays.nnz
Good point.
| inmatrix!(s) || ErrorException("Matrix input.") | ||
| factor!(s) || ErrorException("Numerical Factorization.") | ||
| triangularsolve!(s,rhs) || ErrorException("Triangular Solve.") | ||
| findorder!(s) || error("Finding Order.") |
There was a problem hiding this comment.
Good catch with the exceptions!
|
Does this represent a substantial performance gain? If I recall, symbolic factorization was relatively inexpensive... |
|
This PR is just a robustness fix to the existing implementation of Yup, as far as I know the performance gains from reusing symbolic factorisation are small (10 - 20% from eyeballing profview output?). But this method is already used in eg LinearSolve.jl, so at least from my user's perspective, I think if it exists at all it has to be bulletproof wrt inadvertent changes in sparsity pattern. |
Most definitely. |
|
Is this what you meant to squash: 6c84869 ? |
|
Yup, and I realise there was another wrong turn in the So maybe clearer to just squash-and-merge the whole PR? Or happy to make a new PR with a clean version if you prefer? |
|
New PR might be cleaner, if you don't mind terribly? |
1 similar comment
|
New PR might be cleaner, if you don't mind terribly? |
|
OK done. Closed in favour of cleaned up version in PR #28 Sorry about the mess! |
(update: closed in favour of cleaned up version in PR #28)
Robustness fix to the existing implementation of
sparsepaklu!, adding a check for a change in the sparsity pattern, and addingallow_pattern_changeargument to optionally error if sparsity pattern has changed.Updated docstring now reads:
Calculate LU factorization of a sparse matrix
m, reusing ordering and symbolic factorizationlu, if that was previously calculated.If
allow_pattern_change = true(the default) the sparse matrixmmay have a nonzero pattern different to that of the matrix used to create the LU factorizationlu, in which case the ordering and symbolic factorization are updated.If
allow_pattern_change = falsean error is thrown if the nonzero pattern ofmis different to that of the matrix used to create the LU factorizationlu.If
luhas not been factorized (ie it has just been created with optionfactorize = false) thenluis always updated frommandallow_pattern_changeis ignored.NB: this check is efficient as possible I think, but does need to compare
colptrandrowvalarrays, however this should be inexpensive compared to the actual factorization.Also some bug fixes:
replacing ErrorException() with error()
NB: there is a wrong turn in the first commit (attempting to compare sparsity pattern against the stored graph, which doesn't work as the graph may be modified with an asymmetric matrix), so suggest squash when merging...