Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Deprecate scale #15258

Merged
merged 1 commit into from
Mar 1, 2016
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -128,6 +128,8 @@ Deprecated or removed

* `issym` is deprecated in favor of `issymmetric` to match similar functions (`ishermitian`, ...) ([#15192])

* `scale` is deprecated in favor of either `α*A`, `Diagonal(x)*A`, or `A*Diagonal(x)`. ([#15258])

Julia v0.4.0 Release Notes
==========================

Expand Down
6 changes: 6 additions & 0 deletions base/deprecated.jl
Original file line number Diff line number Diff line change
Expand Up @@ -993,3 +993,9 @@ export call

# Changed issym to issymmetric. #15192
@deprecate issym issymmetric

# 15258
@deprecate scale(α::Number, A::AbstractArray) α*A
@deprecate scale(A::AbstractArray, α::Number) A*α
@deprecate scale(A::AbstractMatrix, x::AbstractVector) A*Diagonal(x)
@deprecate scale(x::AbstractVector, A::AbstractMatrix) Diagonal(x)*A
21 changes: 3 additions & 18 deletions base/docs/helpdb/Base.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4382,20 +4382,6 @@ i-th dimension of `A` should be repeated.
"""
repeat

"""
scale(A, b)
scale(b, A)

Scale an array `A` by a scalar `b`, returning a new array.

If `A` is a matrix and `b` is a vector, then `scale(A,b)` scales each column `i` of `A` by
`b[i]` (similar to `A*diagm(b)`), while `scale(b,A)` scales each row `i` of `A` by `b[i]`
(similar to `diagm(b)*A`), returning a new array.

Note: for large `A`, `scale` can be much faster than `A .* b` or `b .* A`, due to the use of BLAS.
"""
scale

"""
ReentrantLock()

Expand Down Expand Up @@ -6680,12 +6666,11 @@ issetuid
scale!(A, b)
scale!(b, A)

Scale an array `A` by a scalar `b`, similar to [`scale`](:func:`scale`) but overwriting `A`
in-place.
Scale an array `A` by a scalar `b` overwriting `A` in-place.

If `A` is a matrix and `b` is a vector, then `scale!(A,b)` scales each column `i` of `A` by
`b[i]` (similar to `A*diagm(b)`), while `scale!(b,A)` scales each row `i` of `A` by `b[i]`
(similar to `diagm(b)*A`), again operating in-place on `A`.
`b[i]` (similar to `A*Diagonal(b)`), while `scale!(b,A)` scales each row `i` of `A` by `b[i]`
(similar to `Diagonal(b)*A`), again operating in-place on `A`.

"""
scale!
Expand Down
1 change: 0 additions & 1 deletion base/linalg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -111,7 +111,6 @@ export
lqfact!,
lqfact,
rank,
scale,
scale!,
schur,
schurfact!,
Expand Down
4 changes: 2 additions & 2 deletions base/linalg/dense.jl
Original file line number Diff line number Diff line change
Expand Up @@ -179,7 +179,7 @@ function ^(A::Matrix, p::Number)
v, X = eig(A)
any(v.<0) && (v = complex(v))
Xinv = ishermitian(A) ? X' : inv(X)
scale(X, v.^p)*Xinv
(X * Diagonal(v.^p)) * Xinv
end

# Matrix exponential
Expand Down Expand Up @@ -483,7 +483,7 @@ function pinv{T}(A::StridedMatrix{T}, tol::Real)
index = SVD.S .> tol*maximum(SVD.S)
Sinv[index] = one(Stype) ./ SVD.S[index]
Sinv[find(!isfinite(Sinv))] = zero(Stype)
return SVD.Vt'scale(Sinv, SVD.U')
return SVD.Vt' * (Diagonal(Sinv) * SVD.U')
end
function pinv{T}(A::StridedMatrix{T})
tol = eps(real(float(one(T))))*maximum(size(A))
Expand Down
10 changes: 6 additions & 4 deletions base/linalg/diagonal.jl
Original file line number Diff line number Diff line change
Expand Up @@ -106,8 +106,10 @@ end
/{T<:Number}(D::Diagonal, x::T) = Diagonal(D.diag / x)
*(Da::Diagonal, Db::Diagonal) = Diagonal(Da.diag .* Db.diag)
*(D::Diagonal, V::AbstractVector) = D.diag .* V
*(A::AbstractMatrix, D::Diagonal) = scale(A,D.diag)
*(D::Diagonal, A::AbstractMatrix) = scale(D.diag,A)
*(A::AbstractMatrix, D::Diagonal) =
scale!(similar(A, promote_op(MulFun(), eltype(A), eltype(D.diag))), A, D.diag)
*(D::Diagonal, A::AbstractMatrix) =
scale!(similar(A, promote_op(MulFun(), eltype(A), eltype(D.diag))), D.diag, A)

A_mul_B!(A::Diagonal,B::AbstractMatrix) = scale!(A.diag,B)
At_mul_B!(A::Diagonal,B::AbstractMatrix)= scale!(A.diag,B)
Expand Down Expand Up @@ -181,8 +183,8 @@ function A_ldiv_B!(D::Diagonal, B::StridedVecOrMat)
end
return B
end
(\)(D::Diagonal, B::AbstractMatrix) = scale(1 ./ D.diag, B)
(\)(D::Diagonal, b::AbstractVector) = reshape(scale(1 ./ D.diag, reshape(b, length(b), 1)), length(b))
(\)(D::Diagonal, A::AbstractMatrix) = D.diag .\ A
(\)(D::Diagonal, b::AbstractVector) = D.diag .\ b
(\)(Da::Diagonal, Db::Diagonal) = Diagonal(Db.diag ./ Da.diag)

function inv{T}(D::Diagonal{T})
Expand Down
3 changes: 0 additions & 3 deletions base/linalg/generic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,9 +2,6 @@

## linalg.jl: Some generic Linear Algebra definitions

scale(X::AbstractArray, s::Number) = X*s
scale(s::Number, X::AbstractArray) = s*X

# For better performance when input and output are the same array
# See https://github.com/JuliaLang/julia/issues/8415#issuecomment-56608729
function generic_scale!(X::AbstractArray, s::Number)
Expand Down
2 changes: 0 additions & 2 deletions base/linalg/matmul.jl
Original file line number Diff line number Diff line change
Expand Up @@ -38,8 +38,6 @@ function scale!(C::AbstractMatrix, b::AbstractVector, A::AbstractMatrix)
end
C
end
scale(A::AbstractMatrix, b::AbstractVector) = scale!(similar(A, promote_op(MulFun(),eltype(A),eltype(b))), A, b)
scale(b::AbstractVector, A::AbstractMatrix) = scale!(similar(b, promote_op(MulFun(),eltype(b),eltype(A)), size(A)), b, A)

# Dot products

Expand Down
2 changes: 1 addition & 1 deletion base/sparse.jl
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ import Base: @get!, acos, acosd, acot, acotd, acsch, asech, asin, asind, asinh,
exp, expm1, factorize, find, findmax, findmin, findnz, float, full, getindex,
hcat, hvcat, imag, indmax, ishermitian, kron, length, log, log1p, max, min,
maximum, minimum, norm, one, promote_eltype, real, reinterpret, reshape, rot180,
rotl90, rotr90, round, scale, scale!, setindex!, similar, size, transpose, tril,
rotl90, rotr90, round, scale!, setindex!, similar, size, transpose, tril,
triu, vcat, vec

import Base.Broadcast: eltype_plus, broadcast_shape
Expand Down
2 changes: 1 addition & 1 deletion base/sparse/cholmod.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1020,7 +1020,7 @@ function sparse(F::Factor)
else
LD = sparse(F[:LD])
L, d = getLd!(LD)
A = scale(L, d)*L'
A = (L * Diagonal(d)) * L'
end
SparseArrays.sortSparseMatrixCSC!(A)
p = get_perm(F)
Expand Down
6 changes: 0 additions & 6 deletions base/sparse/linalg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -826,12 +826,6 @@ end
scale!(A::SparseMatrixCSC, b::Number) = (scale!(A.nzval, b); A)
scale!(b::Number, A::SparseMatrixCSC) = (scale!(b, A.nzval); A)

scale{Tv,Ti,T}(A::SparseMatrixCSC{Tv,Ti}, b::Vector{T}) =
scale!(similar(A, promote_type(Tv,T)), A, b)

scale{T,Tv,Ti}(b::Vector{T}, A::SparseMatrixCSC{Tv,Ti}) =
scale!(similar(A, promote_type(Tv,T)), b, A)

function factorize(A::SparseMatrixCSC)
m, n = size(A)
if m == n
Expand Down
17 changes: 5 additions & 12 deletions base/sparse/sparsevector.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1191,18 +1191,11 @@ scale!(x::AbstractSparseVector, a::Complex) = (scale!(nonzeros(x), a); x)
scale!(a::Real, x::AbstractSparseVector) = scale!(nonzeros(x), a)
scale!(a::Complex, x::AbstractSparseVector) = scale!(nonzeros(x), a)

scale(x::AbstractSparseVector, a::Real) =
SparseVector(length(x), copy(nonzeroinds(x)), scale(nonzeros(x), a))
scale(x::AbstractSparseVector, a::Complex) =
SparseVector(length(x), copy(nonzeroinds(x)), scale(nonzeros(x), a))

scale(a::Real, x::AbstractSparseVector) = scale(x, a)
scale(a::Complex, x::AbstractSparseVector) = scale(x, a)

*(x::AbstractSparseVector, a::Number) = scale(x, a)
*(a::Number, x::AbstractSparseVector) = scale(x, a)
.*(x::AbstractSparseVector, a::Number) = scale(x, a)
.*(a::Number, x::AbstractSparseVector) = scale(x, a)
# *(x::AbstractSparseVector, a::Real) = SparseVector(length(x), copy(nonzeroinds(x)), scale(nonzeros(x), a))
# *(x::AbstractSparseVector, a::Complex) = SparseVector(length(x), copy(nonzeroinds(x)), scale(nonzeros(x), a))
# *(a::Number, x::AbstractSparseVector) = xa)
# .*(x::AbstractSparseVector, a::Number) = scale(x, a)
# .*(a::Number, x::AbstractSparseVector) = scale(x, a)


# dot
Expand Down
15 changes: 2 additions & 13 deletions doc/stdlib/linalg.rst
Original file line number Diff line number Diff line change
Expand Up @@ -823,25 +823,14 @@ Linear algebra functions in Julia are largely implemented by calling functions f

Construct a diagonal matrix and place ``v`` on the ``k``\ th diagonal.

.. function:: scale(A, b)
scale(b, A)

.. Docstring generated from Julia source

Scale an array ``A`` by a scalar ``b``\ , returning a new array.

If ``A`` is a matrix and ``b`` is a vector, then ``scale(A,b)`` scales each column ``i`` of ``A`` by ``b[i]`` (similar to ``A*diagm(b)``\ ), while ``scale(b,A)`` scales each row ``i`` of ``A`` by ``b[i]`` (similar to ``diagm(b)*A``\ ), returning a new array.

Note: for large ``A``\ , ``scale`` can be much faster than ``A .* b`` or ``b .* A``\ , due to the use of BLAS.

.. function:: scale!(A, b)
scale!(b, A)

.. Docstring generated from Julia source

Scale an array ``A`` by a scalar ``b``\ , similar to :func:`scale` but overwriting ``A`` in-place.
Scale an array ``A`` by a scalar ``b`` overwriting ``A`` in-place.

If ``A`` is a matrix and ``b`` is a vector, then ``scale!(A,b)`` scales each column ``i`` of ``A`` by ``b[i]`` (similar to ``A*diagm(b)``\ ), while ``scale!(b,A)`` scales each row ``i`` of ``A`` by ``b[i]`` (similar to ``diagm(b)*A``\ ), again operating in-place on ``A``\ .
If ``A`` is a matrix and ``b`` is a vector, then ``scale!(A,b)`` scales each column ``i`` of ``A`` by ``b[i]`` (similar to ``A*Diagonal(b)``\ ), while ``scale!(b,A)`` scales each row ``i`` of ``A`` by ``b[i]`` (similar to ``Diagonal(b)*A``\ ), again operating in-place on ``A``\ .

.. function:: Tridiagonal(dl, d, du)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

so is this not true any more?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The deprecation suggestion is to use Diagonal(b)*A which calls BLAS. Actually, the difference is not really because of BLAS, but broadcasting overhead. With a nested loop, the code vectorizes and is as fast as BLAS on my machine.

We could move the "Note" to an entry for *(Diagonal,Matrix) or broadcast, but I'm not sure it will be easy to find in any of those places.


Expand Down
2 changes: 1 addition & 1 deletion test/blas.jl
Original file line number Diff line number Diff line change
Expand Up @@ -113,7 +113,7 @@ for elty in [Float32, Float64, Complex64, Complex128]

# scal
α = rand(elty)
@test BLAS.scal(n,α,a,1) ≈ scale(α,a)
@test BLAS.scal(n,α,a,1) ≈ α * a

# trsv
A = triu(rand(elty,n,n))
Expand Down
4 changes: 2 additions & 2 deletions test/linalg/eigen.jl
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ debug && println("symmetric generalized eigenproblem")
asym_sg = asym[1:n1, 1:n1]
a_sg = a[:,n1+1:n2]
f = eigfact(asym_sg, a_sg'a_sg)
@test_approx_eq asym_sg*f[:vectors] scale(a_sg'a_sg*f[:vectors], f[:values])
@test_approx_eq asym_sg*f[:vectors] (a_sg'a_sg*f[:vectors]) * Diagonal(f[:values])
@test_approx_eq f[:values] eigvals(asym_sg, a_sg'a_sg)
@test_approx_eq_eps prod(f[:values]) prod(eigvals(asym_sg/(a_sg'a_sg))) 200ε
@test eigvecs(asym_sg, a_sg'a_sg) == f[:vectors]
Expand All @@ -66,7 +66,7 @@ debug && println("Non-symmetric generalized eigenproblem")
a1_nsg = a[1:n1, 1:n1]
a2_nsg = a[n1+1:n2, n1+1:n2]
f = eigfact(a1_nsg, a2_nsg)
@test_approx_eq a1_nsg*f[:vectors] scale(a2_nsg*f[:vectors], f[:values])
@test_approx_eq a1_nsg*f[:vectors] (a2_nsg*f[:vectors]) * Diagonal(f[:values])
@test_approx_eq f[:values] eigvals(a1_nsg, a2_nsg)
@test_approx_eq_eps prod(f[:values]) prod(eigvals(a1_nsg/a2_nsg)) 50000ε
@test eigvecs(a1_nsg, a2_nsg) == f[:vectors]
Expand Down
38 changes: 9 additions & 29 deletions test/linalg/generic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -105,26 +105,6 @@ let aa = reshape([1.:6;], (2,3))
a = sub(aa, 1:2, 1:2)
end

# 2-argument version of scale
@test scale(a, 5.) == a*5
@test scale(5., a) == a*5
@test scale([1.; 2.], a) == a.*[1; 2]
@test scale([1; 2], a) == a.*[1; 2]
@test scale(eye(Int, 2), 0.5) == 0.5*eye(2)
@test scale([1; 2], sub(a, :, :)) == a.*[1; 2]
@test scale(sub([1; 2], :), a) == a.*[1; 2]
@test_throws DimensionMismatch scale(ones(3), a)

if atype == "Array"
@test scale(a, [1.; 2.; 3.]) == a.*[1 2 3]
@test scale(a, [1; 2; 3]) == a.*[1 2 3]
@test_throws DimensionMismatch scale(a, ones(2))
else
@test scale(a, [1.; 2.]) == a.*[1 2]
@test scale(a, [1; 2]) == a.*[1 2]
@test_throws DimensionMismatch scale(a, ones(3))
end

# 2-argument version of scale!
@test scale!(copy(a), 5.) == a*5
@test scale!(5., copy(a)) == a*5
Expand Down Expand Up @@ -168,15 +148,15 @@ end

# scale real matrix by complex type
@test_throws InexactError scale!([1.0], 2.0im)
@test isequal(scale([1.0], 2.0im), Complex{Float64}[2.0im])
@test isequal(scale(2.0im, [1.0]), Complex{Float64}[2.0im])
@test isequal(scale(Float32[1.0], 2.0f0im), Complex{Float32}[2.0im])
@test isequal(scale(Float32[1.0], 2.0im), Complex{Float64}[2.0im])
@test isequal(scale(Float64[1.0], 2.0f0im), Complex{Float64}[2.0im])
@test isequal(scale(Float32[1.0], big(2.0)im), Complex{BigFloat}[2.0im])
@test isequal(scale(Float64[1.0], big(2.0)im), Complex{BigFloat}[2.0im])
@test isequal(scale(BigFloat[1.0], 2.0im), Complex{BigFloat}[2.0im])
@test isequal(scale(BigFloat[1.0], 2.0f0im), Complex{BigFloat}[2.0im])
@test isequal([1.0] * 2.0im, Complex{Float64}[2.0im])
@test isequal(2.0im * [1.0], Complex{Float64}[2.0im])
@test isequal(Float32[1.0] * 2.0f0im, Complex{Float32}[2.0im])
@test isequal(Float32[1.0] * 2.0im, Complex{Float64}[2.0im])
@test isequal(Float64[1.0] * 2.0f0im, Complex{Float64}[2.0im])
@test isequal(Float32[1.0] * big(2.0)im, Complex{BigFloat}[2.0im])
@test isequal(Float64[1.0] * big(2.0)im, Complex{BigFloat}[2.0im])
@test isequal(BigFloat[1.0] * 2.0im, Complex{BigFloat}[2.0im])
@test isequal(BigFloat[1.0] * 2.0f0im, Complex{BigFloat}[2.0im])

# test scale and scale! for non-commutative multiplication
q = Quaternion([0.44567, 0.755871, 0.882548, 0.423612])
Expand Down
2 changes: 1 addition & 1 deletion test/linalg/lapack.jl
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ let # syevr
A = convert(Array{elty, 2}, A)
Asym = A'A
vals, Z = LAPACK.syevr!('V', copy(Asym))
@test_approx_eq Z*scale(vals, Z') Asym
@test_approx_eq Z * (Diagonal(vals) * Z') Asym
@test all(vals .> 0.0)
@test_approx_eq LAPACK.syevr!('N','V','U',copy(Asym),0.0,1.0,4,5,-1.0)[1] vals[vals .< 1.0]
@test_approx_eq LAPACK.syevr!('N','I','U',copy(Asym),0.0,1.0,4,5,-1.0)[1] vals[4:5]
Expand Down
2 changes: 1 addition & 1 deletion test/linalg/svd.jl
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ debug && println("\ntype of a: ", eltya, "\n")
debug && println("singular value decomposition")
usv = svdfact(a)
@test usv[:S] === svdvals(usv)
@test usv[:U]*scale(usv[:S],usv[:Vt]) ≈ a
@test usv[:U] * (Diagonal(usv[:S]) * usv[:Vt]) ≈ a
@test full(usv) ≈ a
@test usv[:Vt]' ≈ usv[:V]
@test_throws KeyError usv[:Z]
Expand Down
6 changes: 0 additions & 6 deletions test/linalg/triangular.jl
Original file line number Diff line number Diff line change
Expand Up @@ -207,12 +207,6 @@ for elty1 in (Float32, Float64, BigFloat, Complex64, Complex128, Complex{BigFloa
end
end

@test scale(A1,0.5) == 0.5*A1
@test scale(0.5,A1) == 0.5*A1
@test scale(A1,0.5im) == 0.5im*A1
@test scale(0.5im,A1) == 0.5im*A1


# Binary operations
@test A1*0.5 == full(A1)*0.5
@test 0.5*A1 == 0.5*full(A1)
Expand Down
26 changes: 13 additions & 13 deletions test/sparsedir/sparse.jl
Original file line number Diff line number Diff line change
Expand Up @@ -228,20 +228,20 @@ sA = sprandn(3, 7, 0.5)
sC = similar(sA)
dA = full(sA)
b = randn(7)
@test scale(dA, b) == scale(sA, b)
@test scale(dA, b) == scale!(sC, sA, b)
@test scale(dA, b) == scale!(copy(sA), b)
@test dA * Diagonal(b) == sA * Diagonal(b)
@test dA * Diagonal(b) == scale!(sC, sA, b)
@test dA * Diagonal(b) == scale!(copy(sA), b)
b = randn(3)
@test scale(b, dA) == scale(b, sA)
@test scale(b, dA) == scale!(sC, b, sA)
@test scale(b, dA) == scale!(b, copy(sA))

@test scale(dA, 0.5) == scale(sA, 0.5)
@test scale(dA, 0.5) == scale!(sC, sA, 0.5)
@test scale(dA, 0.5) == scale!(copy(sA), 0.5)
@test scale(0.5, dA) == scale(0.5, sA)
@test scale(0.5, dA) == scale!(sC, sA, 0.5)
@test scale(0.5, dA) == scale!(0.5, copy(sA))
@test Diagonal(b) * dA == Diagonal(b) * sA
@test Diagonal(b) * dA == scale!(sC, b, sA)
@test Diagonal(b) * dA == scale!(b, copy(sA))

@test dA * 0.5 == sA * 0.5
@test dA * 0.5 == scale!(sC, sA, 0.5)
@test dA * 0.5 == scale!(copy(sA), 0.5)
@test 0.5 * dA == 0.5 * sA
@test 0.5 * dA == scale!(sC, sA, 0.5)
@test 0.5 * dA == scale!(0.5, copy(sA))
@test scale!(sC, 0.5, sA) == scale!(sC, sA, 0.5)

# copy!
Expand Down
8 changes: 4 additions & 4 deletions test/sparsedir/sparsevector.jl
Original file line number Diff line number Diff line change
Expand Up @@ -580,10 +580,10 @@ let x = sprand(16, 0.5), x2 = sprand(16, 0.4)

# scale
let sx = SparseVector(x.n, x.nzind, x.nzval * 2.5)
@test exact_equal(scale(x, 2.5), sx)
@test exact_equal(scale(x, 2.5 + 0.0*im), complex(sx))
@test exact_equal(scale(2.5, x), sx)
@test exact_equal(scale(2.5 + 0.0*im, x), complex(sx))
@test exact_equal(x * 2.5, sx)
@test exact_equal(x * (2.5 + 0.0*im), complex(sx))
@test exact_equal(2.5 * x, sx)
@test exact_equal((2.5 + 0.0*im) * x, complex(sx))
@test exact_equal(x * 2.5, sx)
@test exact_equal(2.5 * x, sx)
@test exact_equal(x .* 2.5, sx)
Expand Down
6 changes: 3 additions & 3 deletions test/sparsedir/umfpack.jl
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ for Tv in (Float64, Complex128)
lua = lufact(A)
@test nnz(lua) == 18
L,U,p,q,Rs = lua[:(:)]
@test_approx_eq scale(Rs,A)[p,q] L*U
@test_approx_eq (Diagonal(Rs) * A)[p,q] L * U

@test_approx_eq det(lua) det(full(A))

Expand All @@ -45,15 +45,15 @@ for Ti in Base.SparseArrays.UMFPACK.UMFITypes.types
Ac = convert(SparseMatrixCSC{Complex128,Ti}, Ac0)
lua = lufact(Ac)
L,U,p,q,Rs = lua[:(:)]
@test_approx_eq scale(Rs,Ac)[p,q] L*U
@test_approx_eq (Diagonal(Rs) * Ac)[p,q] L * U
end

for elty in (Float64, Complex128)
for (m, n) in ((10,5), (5, 10))
A = sparse([1:min(m,n); rand(1:m, 10)], [1:min(m,n); rand(1:n, 10)], elty == Float64 ? randn(min(m, n) + 10) : complex(randn(min(m, n) + 10), randn(min(m, n) + 10)))
F = lufact(A)
L, U, p, q, Rs = F[:(:)]
@test_approx_eq scale(Rs,A)[p,q] L*U
@test_approx_eq (Diagonal(Rs) * A)[p,q] L * U
end
end

Expand Down