In-place operations are incorrect on UnitTriangular inputs

[tkelman]

Due to calling triu!(B.data) or tril!(B.data) without setting B's diagonals to 1 first. Ref #16205 (comment)

 _/ |\__'_|_|_|\__'_|  |  Commit c2d9502* (0 days old master)
|__/                   |  x86_64-w64-mingw32

julia> A = UpperTriangular(rand(4,4))
4×4 UpperTriangular{Float64,Array{Float64,2}}:
 0.667194  0.707125  0.45711   0.273548
  ⋅        0.13787   0.307034  0.600464
  ⋅         ⋅        0.314295  0.113895
  ⋅         ⋅         ⋅        0.719889

julia> B = LinAlg.UnitUpperTriangular(rand(4,4))
4×4 Base.LinAlg.UnitUpperTriangular{Float64,Array{Float64,2}}:
 1.0  0.877432  0.691445  0.100324
 0.0  1.0       0.405473  0.568657
 0.0  0.0       1.0       0.634321
 0.0  0.0       0.0       1.0

julia> B.data
4×4 Array{Float64,2}:
 0.497388  0.877432  0.691445  0.100324
 0.615141  0.645568  0.405473  0.568657
 0.18405   0.206572  0.164716  0.634321
 0.118562  0.886125  0.335527  0.0684593

julia> full(A) * full(B)
4×4 Array{Float64,2}:
 0.667194  1.29254  1.20516   1.03255
 0.0       0.13787  0.362937  0.873623
 0.0       0.0      0.314295  0.313259
 0.0       0.0      0.0       0.719889

julia> C = LinAlg.A_mul_B!(A, B)
4×4 UpperTriangular{Float64,Array{Float64,2}}:
 0.331854  1.04192    0.823342   0.777728
  ⋅        0.0890047  0.106476   0.314266
  ⋅         ⋅         0.0517696  0.207161
  ⋅         ⋅          ⋅         0.0492832

[tkelman]

cc @andreasnoack