Dividing diagonal by triangular matrix returns sparse matrix

[hbarthels]

Depending on what types are used, dividing a diagonal matrix by a triangular matrix returns a sparse matrix, even though the result is triangular. I think it would make a lot more sense to return a matrix with a triangular type instead.

julia> VERSION
v"1.3.0-rc1.0"

julia> n = 3;

julia> L = rand(n, n);

julia> D = rand(n);

julia> tril(L)\diagm(D)
3×3 Array{Float64,2}:
  0.0666276   0.0       0.0     
 -0.115247    1.0267    0.0     
  0.0406685  -0.656616  0.507579

julia> LowerTriangular(L)\diagm(D)
3×3 Array{Float64,2}:
  0.0666276   0.0       0.0     
 -0.115247    1.0267    0.0     
  0.0406685  -0.656616  0.507579

julia> tril(L)\Diagonal(D)
3×3 SparseArrays.SparseMatrixCSC{Float64,Int64} with 6 stored entries:
  [1, 1]  =  0.0666276
  [2, 1]  =  -0.115247
  [3, 1]  =  0.0406685
  [2, 2]  =  1.0267
  [3, 2]  =  -0.656616
  [3, 3]  =  0.507579

julia> LowerTriangular(L)\Diagonal(D)
3×3 SparseArrays.SparseMatrixCSC{Float64,Int64} with 6 stored entries:
  [1, 1]  =  0.0666276
  [2, 1]  =  -0.115247
  [3, 1]  =  0.0406685
  [2, 2]  =  1.0267
  [3, 2]  =  -0.656616
  [3, 3]  =  0.507579

It's very similar to what was addressed here: #27999

[brenhinkeller]

Looks like all of these now (1.8.3) return dense matrices

julia> tril(L)\diagm(D)
3×3 Matrix{Float64}:
  0.63596    0.0       0.0
 -0.294611   0.524089  0.0
  0.193844  -1.02639   0.375219

julia> LowerTriangular(L)\diagm(D)
3×3 Matrix{Float64}:
  0.63596    0.0       0.0
 -0.294611   0.524089  0.0
  0.193844  -1.02639   0.375219

julia> tril(L)\Diagonal(D)
3×3 Matrix{Float64}:
  0.63596    0.0       0.0
 -0.294611   0.524089  0.0
  0.193844  -1.02639   0.375219

julia> LowerTriangular(L)\Diagonal(D)
3×3 Matrix{Float64}:
  0.63596    0.0       0.0
 -0.294611   0.524089  0.0
  0.193844  -1.02639   0.375219

There still may be a performance improvement on the table here though if we could return Triangular rather than full dense matrices