Skip to content

Conversation

@dpo
Copy link
Contributor

@dpo dpo commented Apr 26, 2020

This PR adds the ordering kwarg to qr() for sparse matrices (which calls SuiteSparseQR) so the user can select one of the orderings exposed by SPQR.
The default value of ordering is ORDERING_DEFAULT so that the default behavior is unchanged.

A little care must be taken when returning from _ldiv_basic because if ordering == ORDERING_FIXED, SPQR returns an empty column permutation (to indicate that no permutation occurred). I tried writing the return statement so it's more readable but got errors in the threaded SuiteSparse tests.

I added tests for over- and under-determined problems.

@andreasnoack

@dkarrasch dkarrasch added linear algebra Linear algebra sparse Sparse arrays labels Apr 26, 2020
# Apply right permutation and extract solution from X
return getindex(X, ntuple(i -> i == 1 ? invperm(F.cpiv) : :, Val(ndims(B)))...)
# NB: cpiv == [] if SPQR was called with ORDERING_FIXED
return length(F.cpiv) == 0 ? getindex(X, ntuple(i -> i == 1 ? (1:size(F,2)) : :, Val(ndims(B)))...) : getindex(X, ntuple(i -> i == 1 ? invperm(F.cpiv) : :, Val(ndims(B)))...)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Style nit: please keep line length reasonable. So usually about 1 anonymous function per expression.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I had trouble with this when I tested. That's what I try to explain above. Adding a return statement for the special case causes threaded tests to fail. I'm not sure why. I'll try again.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh well, nevermind. I cleaned it up and it seems ok.

@ViralBShah
Copy link
Member

Specifying the ordering as a value from an array is a bit difficult to use. Perhaps use :Ordering_name. We should add all the orderings to the documentation, and also a NEWS entry.

@dpo
Copy link
Contributor Author

dpo commented Apr 28, 2020

You don't refer to them as a value in an array. You refer to them by name: ORDERING_AMD, ORDERING_METIS, etc.

@ViralBShah
Copy link
Member

Ah thanks. I didn't notice that was the case only in ErrorException.

@andreasnoack andreasnoack merged commit c8797ad into JuliaLang:master Apr 30, 2020
@dpo
Copy link
Contributor Author

dpo commented Apr 30, 2020

Sorry I forgot to update docs and news. Will do in a separate PR.

@dpo
Copy link
Contributor Author

dpo commented May 1, 2020

Could you indicate where I should modify the docs? Every time I read https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.qr I'm surprised that it doesn't seem to talk about sparse QR at all, but I can't find another place that does. The sparse QR doesn't take a blocksize argument, and the factorization object doesn't have the p and P attributes (it has pcol and prow).

@ViralBShah
Copy link
Member

Yes, unfortunately none of the sparse factorizations have documentation. I think the right place to put them is in the sparse matrix implementation, and the doc system should find it. At least that's my understanding.

@fredrikekre should be able to confirm that.

@andreasnoack
Copy link
Member

All the sparse factorizations have docstrings but it does look like that they aren't included in the documentation. I guess we need to generate a file similar to https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/docs/src/index.md in the stdlib/SuiteSparse folder.

@ViralBShah
Copy link
Member

@andreasnoack Can you open a separate issue for that?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

linear algebra Linear algebra sparse Sparse arrays

Projects

None yet

Development

Successfully merging this pull request may close these issues.

5 participants