Preserve structured matrix types more often when broadcasting

[mbauman]

The primary motivation here is exponentiation (fixes #32759). Previously X .^ 2 would return a Matrix but two=2; X .^ two would return the same structured matrix type as X. This simply teaches LinearAlgebra a little bit more about broadcast so it can safely compute the zero-preservation property of broadcasts involving Refs, which is sufficient to handle the literal pow argument list.

[mbauman]

Let's re-run CI and get this in.

[mbauman]

Interesting — the 32 bit failures look related. I don't see any differences in the printout below, but I think the � is a mojibake ⋯.

Test Failed at C:\cygwin\home\Administrator\buildbot\worker-tabularasa\tester_win32\build\share\julia\stdlib\v1.3\LinearAlgebra\test\structuredbroadcast.jl:37
  Expression: X .^ 2 == X .^ (2,) == fX .^ 2 == X .^ two
   Evaluated: 
    [0.31502726237096007 0.43928288933667414 � 0.0 0.0; 0.41656112475245854 0.009615577025932596 � 0.0 0.0; � ; 0.0 0.0 � 0.037316269381872484 0.2966281365862301; 0.0 0.0 � 0.788465378213704 0.7921566497473703] ==
    [0.31502726237096007 0.43928288933667414 � 0.0 0.0; 0.41656112475245854 0.009615577025932596 � 0.0 0.0; � ; 0.0 0.0 � 0.037316269381872484 0.2966281365862301; 0.0 0.0 � 0.788465378213704 0.7921566497473703] == 
    [0.31502726237096007 0.43928288933667414 � 0.0 0.0; 0.41656112475245854 0.009615577025932596 � 0.0 0.0; � ; 0.0 0.0 � 0.037316269381872484 0.2966281365862301; 0.0 0.0 � 0.788465378213704 0.7921566497473703] == 
    [0.31502726237096007 0.43928288933667414 � 0.0 0.0; 0.41656112475245854 0.009615577025932596 � 0.0 0.0; � ; 0.0 0.0 � 0.037316269381872484 0.2966281365862301; 0.0 0.0 � 0.788465378213704 0.7921566497473703]

[dkarrasch]

@mbauman Could you rebase this PR? The last CI failures look weird, and partly unrelated. Would be nice to see if those issues have been resolved meanwhile, so we could get this in.

[mbauman]

Bah, now it passes tests? Looks like this is a failure that depends upon the RNG seed — likely due to differences between BLAS' exponentiation and the structured optimization. We're probably within an ULP or two. I suppose I'll just make those ≈s.

[mbauman]

Wait that makes no sense. We're not hitting BLAS in any case here. Not sure where a difference would come in.

[dkarrasch]

Could it be that it has to do with very small numbers getting squared and getting even smaller? It has to be somehow related to the relative accuracy, I guess. What if you define

D = Diagonal(rand(N).+1)

or

D = Diagonal(rand(-10:10, N))

?

[mbauman]

This is fun:

julia> two = 2
2

julia> x = 0.9390576568834563
0.9390576568834563

julia> x ^ 2
0.8818292829514472

julia> x ^ two
0.8818292829514471

Edit: looks like the difference is between using LLVM's pow intrinsic vs. fmul doubling:

julia> @code_llvm Base.literal_pow(^, x, Val(2))

;  @ intfuncs.jl:244 within `literal_pow'
define double @julia_literal_pow_19302(double) {
top:
; ┌ @ float.jl:405 within `*'
   %1 = fmul double %0, %0
; └
  ret double %1
}

julia> @code_llvm x^2

;  @ math.jl:860 within `^'
define double @"julia_^_19325"(double, i64) {
top:
; ┌ @ float.jl:60 within `Float64'
   %2 = sitofp i64 %1 to double
; └
  %3 = call double @llvm.pow.f64(double %0, double %2)
  ret double %3
}

[dkarrasch]

Which one corresponds to which? x^2 is the first, and x^two the second? Like literal vs. variable?

[dkarrasch]

The failure is in Distributed, I'd say we can merge this, right? @mbauman

[mbauman]

Since #34863 got in first, I think we should drop its workaround commit from here.

[dkarrasch]

Oh yea, you mean the "fine-grained tests" thing?

[mbauman]

Sure, that one too. There had also been another commit that completely avoided looking at numeric equality between X .^ 2 and X .^ two, but now those are the same.