Add missing optimization for *,/ between Diagonal and Triangular

[N5N3]

This PR tries to add missing optimization for *, / related fuctions between Diagonal and AbstractTriangular.
It was first discussed in #42321 (comment), and I notice that such patten also applies to /, ldiv!, \ and rdiv!.
I believe that unified implement help us cover all possible optimization.
For example, lmul! has been missing for years.

Test has been added.

(I hope I didn't ignore anything.)

[N5N3]

Does failure on freebsd64 related to this PR?
I see qr and triangular, but it passed once before I added the require_one_based_indexing(B).
It seems win32 and win64's tests are not running?

[dkarrasch]

freebsd64 can be ignored, if this is the only one failing.

[dkarrasch]

A few more comments. Otherwise, this is shaping up nicely.

[N5N3]

Error found in

@test Diagonal(spzeros(5)) \ view(rand(10), 1:5) == [Inf, Inf, Inf, Inf, Inf]

This throw SingularException(1) after this PR.
Just change it to @test_throws ?

Others are unexpected pass

@test H/x == Array(H)/x broken = eltype(H) <: Furlong && x isa Union{Bidiagonal, Diagonal, UpperTriangular}
@test x\H == x\Array(H) broken = eltype(H) <: Furlong && x isa Union{Bidiagonal, Diagonal, UpperTriangular}
@test H/x isa UpperHessenberg broken = eltype(H) <: Furlong && x isa Union{Bidiagonal, Diagonal}
@test x\H isa UpperHessenberg broken = eltype(H) <: Furlong && x isa Union{Bidiagonal, Diagonal}

I haven't look into the reason, but all these are about Diagonal.
Is broken the excepted behaviour? Or we just remove Diagonal here.

[dkarrasch]

This throw SingularException(1) after this PR.
Just change it to @test_throws ?

I'd say yes, but one could argue that diagonal solves are nothing but a sequence of scalar divisions. When we have scalar 5/0 we get Inf without a SingularException or any other warning. In any case, we need to make sure that

@test Diagonal(spzeros(5)) \ view(rand(10), 1:5) == Matrix(Diagonal(spzeros(5))) \ view(rand(10), 1:5)

If we change this for diagonal, this is going to be passed on the the matrix one, because it checks for triangularity:

function (\)(A::AbstractMatrix, B::AbstractVecOrMat)
    require_one_based_indexing(A, B)
    m, n = size(A)
    if m == n
        if istril(A)
            if istriu(A)
                return Diagonal(A) \ B
...

I'd vote for consistency with

julia> lu(zeros(5,5))
ERROR: SingularException(1)

so with what you have.

As for the second issue, I guess you resolved broken cases along the way. Congrats! 🎉 🤣

[dkarrasch]

Since this PR almost preserves the number of lines of code, I think it is acceptable to go the last bit and avoid the vec(permutedims construct by some straightforward loop based code. I especially like the fact that this removes a couple of broken tests. Really awesome!

[dkarrasch]

Up to one last suggestion, this is ready for review. The only potentially controversial aspect that I see is the change from "blind" division to a checked solve with throwing SingularExceptions.

[N5N3]

The failure on tester_win32 shows that:

      From worker 8:	
      From worker 8:	Please submit a bug report with steps to reproduce this fault, and any error messages that follow (in their entirety). Thanks.
      From worker 8:	Exception: EXCEPTION_ACCESS_VIOLATION at 0x156d8264 -- .text at C:\buildbot\worker-tabularasa\tester_win32\build\bin\libopenblas.dll (unknown line)
      From worker 8:	in expression starting at C:\buildbot\worker-tabularasa\tester_win32\build\share\julia\stdlib\v1.8\LinearAlgebra\test\tridiag.jl:28
      From worker 8:	.text at C:\buildbot\worker-tabularasa\tester_win32\build\bin\libopenblas.dll (unknown line)
      From worker 8:	.text at C:\buildbot\worker-tabularasa\tester_win32\build\bin\libopenblas.dll (unknown line)
      From worker 8:	.text at C:\buildbot\worker-tabularasa\tester_win32\build\bin\libopenblas.dll (unknown line)
      From worker 8:	stegr! at C:\buildbot\worker\package_win32\build\usr\share\julia\stdlib\v1.8\LinearAlgebra\src\lapack.jl:3853
      From worker 8:	eigen! at C:\buildbot\worker\package_win32\build\usr\share\julia\stdlib\v1.8\LinearAlgebra\src\tridiag.jl:294 [inlined]
      From worker 8:	eigen at C:\buildbot\worker\package_win32\build\usr\share\julia\stdlib\v1.8\LinearAlgebra\src\tridiag.jl:296
      From worker 8:	unknown function (ip: 016c7eed)

This seems to be an openblas issue?

[dkarrasch]

I'm not sure I like the _promote_dotop thing, but otherwise this LGTM. ;-)

[N5N3]

I add 5-arg mul! optimization between Triangular and Diagonal.
The problem is how to handle UnitUpper/Lower input.
Currently, the output's diagonal is first cached in to a diag′::Vector if β != 0, so the speed of 3-arg mul! should not be affected.
For 5-arg mul! the influence is small at least for N >= 100.
Some benchmarks:

n = 100
UTA = UnitUpperTriangular(zeros(n,n));
TA = UpperTriangular(zeros(n,n));
D = Diagonal(zeros(n));
B = TA * D;
@btime mul!($B,$D,$UTA);    # 1.640 μs (0 allocations: 0 bytes)
@btime mul!($B,$D,$TA);       #   1.530 μs (0 allocations: 0 bytes)
@btime mul!($B,$D,$UTA, $0, $1);   #   5.867 μs (5 allocations: 1008 bytes)
@btime mul!($B,$D,$TA, $0, $1);    #   5.567 μs (3 allocations: 80 bytes)
@btime mul!($B,$D,$UTA, $2, $1);   #   7.250 μs (6 allocations: 1.03 KiB)
@btime mul!($B,$D,$TA, $2, $1);   #   6.900 μs (4 allocations: 128 bytes)

[dkarrasch]

@andreasnoack Do you think there's anything wrong about this change in behavior? #42343 (comment)

Absent objections, I'd merge this soon.