Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Make *Triangular handle units #43972

Merged
merged 34 commits into from
May 11, 2023
Merged

Make *Triangular handle units #43972

merged 34 commits into from
May 11, 2023

Conversation

dkarrasch
Copy link
Member

@dkarrasch dkarrasch commented Jan 28, 2022

This is an attempt to make triangular linear algebra handle unitful objects well, like (most of) UpperHessenberg and Bidiagonal. The first commit rewrites the generic functions in terms of 3-arg mul!, ldiv! and the internal _rdiv! and doesn't change any tests to convince ourselves that this doesn't change behavior. 2-arg versions for BLAS-able argument types remain fully functional. Some of the loops had to be reorder to keep the variables consistent with units.

Closes #42550.

EDIT (copied from below):

The features of the PR are:

  • multiplication and solving with triangular matrices involving units (the triangular matrix, the other factor, or both) works
  • inversion of unitful triangular matrices works
  • matrix square root of unitful UnitTriangular matrices works
  • right in-place multiplication by a triangular matrix does no longer fall back to generic code (see mul!(::Matrix{T}, ::Matrix{T}, ::AbstractTriangular{T, <:Matrix{T}}) where T<:BlasFloat does not use BLAS #42550)
  • 3-arg triangular multiplication/division of triangular matrices is directed towards BLAS if possible
  • reduces mul! methods down to 55 (we have 138 in v1.9)
  • there is a "_ustrip" mechanism that packages can hook into to have generic substitution codes run without overhead by division by some oneunit

@dkarrasch
Copy link
Member Author

@nanosoldier runbenchmarks("linalg", vs = ":master")

@dkarrasch
Copy link
Member Author

@nanosoldier runtests(ALL, vs = ":master")

@nanosoldier
Copy link
Collaborator

Your benchmark job has completed - possible performance regressions were detected. A full report can be found here.

@nanosoldier
Copy link
Collaborator

Your package evaluation job has completed - possible new issues were detected. A full report can be found here.

end
c[i] = A[i,i] \ bi
Copy link
Contributor

@KlausC KlausC Mar 9, 2022

Choose a reason for hiding this comment

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

Isn't there a more efficient way to move the correct physical dimensions to bi?

When the element types of A and b are unitful integers, the result could also be unitful integers, without any floating point operations performed (no division by 1).

Copy link
Member Author

Choose a reason for hiding this comment

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

I don't know of any. convert and alike fail because you can't convert, for instance, square meters to meters or a unitless number. Also, I'm not sure if there is a generic way of identifying whether a unitful quantity is integer-based. For unitless numbers you can see that I made an effort to preserve the integer return type.

Copy link
Member Author

Choose a reason for hiding this comment

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

I worked around the issue. First, the target c is initialized with an integer unitful eltype, and second I replaced the division by a multiplication. I'll add a test for that also.

Copy link
Member Author

Choose a reason for hiding this comment

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

Actually, no. For unitful quantities, we need to "perform" the division to get the right unit: assume both A and b have integer units of kilometers, then the result is "unit-free". In terms of types, this means (Int*km) \ (Int*km) = Float64, otherwise the kilometers don't cancel out.

@dkarrasch
Copy link
Member Author

@nanosoldier runtests(ALL, vs = ":master")

@dkarrasch
Copy link
Member Author

what do you need to make progress here?

Feedback. 😜

This is a pretty big refactoring, and it increases CI runtime in the triangular test file (due to compilation of more methods) and in the bidiag test file due to quite some more tests involving units. I'm wondering if this is all worth it or if we have other options to handle unitful linear algebra (and whether this is all of sufficient interest).

@dkarrasch
Copy link
Member Author

@nanosoldier runtests()

@nanosoldier
Copy link
Collaborator

Your package evaluation job has completed - possible new issues were detected.
A full report can be found here.

@dkarrasch
Copy link
Member Author

@nanosoldier runtests(["ReusableFunctions", "Pyehtim", "CommutativeRings", "HomotopyOpt", "GradientRobustMultiPhysics", "HTTP", "JuMP", "LazyBandedMatrices", "SignalTablesInterface_CairoMakie", "Pickle", "SignalTablesInterface_WGLMakie", "CompressiveLearning", "SummationByPartsOperators", "UnfoldMakie", "Decapodes", "Estapir", "LiteQTL", "OptimizationMetaheuristics", "Kronecker", "ClimaCore"])

@nanosoldier
Copy link
Collaborator

Your package evaluation job has completed - possible new issues were detected.
A full report can be found here.

@dkarrasch
Copy link
Member Author

pkgeval says we're good. So let's go?

@dkarrasch
Copy link
Member Author

@nanosoldier runbenchmarks()

@dkarrasch
Copy link
Member Author

@nanosoldier runbenchmarks("linalg")

@nanosoldier
Copy link
Collaborator

Your benchmark job has completed - successfully executed benchmarks. A full report can be found here.

@dkarrasch
Copy link
Member Author

@nanosoldier runbenchmarks("linalg")

@nanosoldier
Copy link
Collaborator

Your benchmark job has completed - successfully executed benchmarks. A full report can be found here.

@dkarrasch
Copy link
Member Author

I must have missed some update for nanosoldier. We used to have green checkmarks and red crosses in the benchmarks output. How do I read the output now? Does the absence of these signs mean "everything within expectation"? @maleadt @KristofferC

@maleadt
Copy link
Member

maleadt commented May 11, 2023

That's because you didn't do a comparison benchmark. Only the PkgEval Nanosoldier instance supports leaving off the vs keyword.
@vtjnash This has happened a couple of times now. Can you rebase your Nanosoldier instance on top of the latest Nanosoldier.jl?

@dkarrasch
Copy link
Member Author

I wasn't reading carefully, but just looked at the examples at https://github.com/JuliaCI/Nanosoldier.jl, and copy-pasted. Now I see way below that there are more details to be specified. Let's try again.

@nanosoldier runbenchmarks("linalg", vs = ":master")

@nanosoldier
Copy link
Collaborator

Your benchmark job has completed - possible performance regressions were detected. A full report can be found here.

@dkarrasch
Copy link
Member Author

The only reported regression must be fake. The code path (mul! for 2x2 matrices) is not touched here, and I also can't reproduce locally.

@dkarrasch dkarrasch merged commit b21f100 into master May 11, 2023
@dkarrasch dkarrasch deleted the dk/triangular branch May 11, 2023 19:06
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging this pull request may close these issues.

mul!(::Matrix{T}, ::Matrix{T}, ::AbstractTriangular{T, <:Matrix{T}}) where T<:BlasFloat does not use BLAS
6 participants