Fewer MulAddMul branches in Diagonal-triangular mul (#1272)

This reduces TTFX in `Diagonal` - triangular multiplications.
```julia
julia> using Random, LinearAlgebra

julia> A = rand(4,4); U = UpperTriangular(A); D = Diagonal(A);

julia> @time D * U;
  0.131110 seconds (239.39 k allocations: 12.162 MiB, 99.95% compilation time) # master
  0.102569 seconds (227.44 k allocations: 11.472 MiB, 99.94% compilation time) # this PR
```
If the `Diagonal` is on the right, the TTFX is almost identical, but
allocations go down slightly.
```julia
julia> @time U * D;
  0.125025 seconds (221.82 k allocations: 11.252 MiB, 99.95% compilation time) # master
  0.127002 seconds (215.59 k allocations: 10.938 MiB, 12.06% gc time, 99.95% compilation time) # this PR
```

Author: jishnub

src/diagonal.jl

@@ -470,19 +470,19 @@ function __muldiag_nonzeroalpha!(out, D::Diagonal, B::UpperOrLowerTriangular, al
     for j in axes(B, 2)
         # store the diagonal separately for unit triangular matrices
         if isunit
-            @inbounds @stable_muladdmul _modify!(MulAddMul(alpha,beta), D.diag[j] * B[j,j], out, (j,j))
+            @inbounds _modify_nonzeroalpha!(D.diag[j] * B[j,j], out, (j,j), alpha, beta)
         end
         # The indices of out corresponding to the stored indices of B
         rowrange = _rowrange_tri_stored(B, j)
         @inbounds @simd for i in rowrange
-            @stable_muladdmul _modify!(MulAddMul(alpha,beta), D.diag[i] * B_maybeparent[i,j], out_maybeparent, (i,j))
+            _modify_nonzeroalpha!(D.diag[i] * B_maybeparent[i,j], out_maybeparent, (i,j), alpha, beta)
         end
         # Fill the indices of out corresponding to the zeros of B
         # we only fill these if out and B don't have matching zeros
         if !_has_matching_zeros(out, B)
             rowrange = _rowrange_tri_zeros(B, j)
             @inbounds @simd for i in rowrange
-                @stable_muladdmul _modify!(MulAddMul(alpha,beta), D.diag[i] * B[i,j], out, (i,j))
+                _modify_nonzeroalpha!(D.diag[i] * B[i,j], out, (i,j), alpha, beta)
             end
         end
     end
@@ -511,7 +511,7 @@ function __muldiag_nonzeroalpha!(out, A::UpperOrLowerTriangular, D::Diagonal, al
     # we may directly read and write from the parents
     out_maybeparent, A_maybeparent = _has_matching_zeros(out, A) ? (parent(out), parent(A)) : (out, A)
     for j in axes(A, 2)
-        dja = @stable_muladdmul MulAddMul(alpha,false)(@inbounds D.diag[j])
+        dja = @inbounds _djalpha_nonzero(D.diag[j], alpha)
         # store the diagonal separately for unit triangular matrices
         if isunit
             # since alpha is multiplied to the diagonal element of D,
@@ -547,7 +547,7 @@ end
     d2 = D2.diag
     outd = out.diag
     @inbounds @simd for i in eachindex(d1, d2, outd)
-        @stable_muladdmul _modify!(MulAddMul(alpha,beta), d1[i] * d2[i], outd, i)
+        _modify_nonzeroalpha!(d1[i] * d2[i], outd, i, alpha, beta)
     end
     return out
 end