-
-
Notifications
You must be signed in to change notification settings - Fork 5.5k
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
More tests for Triangular, div and mult #31831
Conversation
Would you mind changing the definitions to long form functions now that you looking into these anyway? They are really hard to read in their current form. |
Can do, and I agree about readability ... this pr took a long time to make 😢 |
@@ -357,11 +374,16 @@ for elty1 in (Float32, Float64, BigFloat, ComplexF32, ComplexF64, Complex{BigFlo | |||
@test transpose(A1)*B ≈ transpose(Matrix(A1))*B | |||
@test A1'B ≈ Matrix(A1)'B | |||
@test A1*transpose(B) ≈ Matrix(A1)*transpose(B) | |||
@test A1*Transpose(B) ≈ Matrix(A1)*transpose(B) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
How is this one different from the one above?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I guess it's not. IIUC then B
is a dense matrix without extra structure, and transpose(B)
calls the wrapper Transpose(B)
.
(B = transB.parent; LowerTriangular(rdiv!(tril!(A.data), transpose(B)))) | ||
function lmul!(adjA::Adjoint{<:Any,<:Union{LowerTriangular,UnitLowerTriangular}}, B::UpperTriangular) | ||
A = adjA.parent | ||
return UpperTriangular(lmul!(adjoint(A), triu!(B.data))) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is one of these cases that leave me breathlessly confused. Under which circumstances could adjoint(adjA.parent)
be different from adjA
in the first place? Why isn't this simply
lmul!(adjA::Adjoint{<:Any,<:Union{LowerTriangular,UnitLowerTriangular}}, B::UpperTriangular) =
UpperTriangular(lmul!(adjA, triu!(B.data)))
?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I suspect it might be a consequence of a rather mechanical conversion from the old Ac_mul_B
to the new Adjoint
based. There were a lot of methods to change so I'm not surprised if we can simplify things along the way.
@@ -329,6 +331,21 @@ for elty1 in (Float32, Float64, BigFloat, ComplexF32, ComplexF64, Complex{BigFlo | |||
@test_throws DimensionMismatch transpose(A2) * offsizeA | |||
@test_throws DimensionMismatch A2' * offsizeA | |||
@test_throws DimensionMismatch A2 * offsizeA | |||
if( uplo1 == uplo2 && elty1 == elty2 != Int && t1 != UnitLowerTriangular && t1 != UnitUpperTriangular ) | |||
@test rdiv!(copy(A1), copy(A2)) ≈ A1/A2 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm not sure, but don't you want to turn the rhs into dense matrices? A1/A2
uses rdiv!
, so you may end up testing the same function on both sides. But it's very hard to tell which methods are used eventually.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
All three, A1
, A2
and B
get wrapped by Transpose/Adjoint
when you call transpose
or adjoint
on them:
julia/stdlib/LinearAlgebra/src/triangular.jl
Lines 376 to 383 in 93c9ae4
adjoint(A::LowerTriangular) = Adjoint(A) | |
adjoint(A::UpperTriangular) = Adjoint(A) | |
adjoint(A::UnitLowerTriangular) = Adjoint(A) | |
adjoint(A::UnitUpperTriangular) = Adjoint(A) | |
transpose(A::LowerTriangular) = Transpose(A) | |
transpose(A::UpperTriangular) = Transpose(A) | |
transpose(A::UnitLowerTriangular) = Transpose(A) | |
transpose(A::UnitUpperTriangular) = Transpose(A) |
So there is no need for a distinction or double testing. So in the end, I think it would suffice to use only adjoint
and transpose
.
f98b7bb
to
40e4473
Compare
@kshyatt I hope you don't mind that I resumed your work here. As per the discussion above, I simplified the multiplication code slightly, removed some unnecessary tests, and made sure we test against independent |
"18 successful checks"... I love it! |
If there are no objections, I'd merge this tonight. |
Sounds good to me! |
I fixed a few more cases where we unwrap and immediately wrap again. There are, on the other hand, many cases where the unwrapping makes sense: after unwrapping, we promote the eltype, convert the parent array, and only after that we hit it with the appropriate |
* More tests for Triangular, div and mult * Long form function defs * avoid useless unwrapping-wrapping in triangular multiplication Co-authored-by: Daniel Karrasch <[email protected]>
* More tests for Triangular, div and mult * Long form function defs * avoid useless unwrapping-wrapping in triangular multiplication Co-authored-by: Daniel Karrasch <[email protected]>
A lot of these aren't covered. There are a lot of conditions on the matrix formats and making sure you don't try to sneakily turn a real matrix complex.