How to use ldiv! correctly?

General Usage
1

How to solve the below error when using ldiv!?

julia> Sn=[1 2 3;0 4 0;0 0 1]
3×3 Matrix{Int64}:
 1  2  3
 0  4  0
 0  0  1

julia> Rn =[1,2,3]
3-element Vector{Int64}:
 1
 2
 3

julia> Qn= zeros(3)
3-element Vector{Float64}:
 0.0
 0.0
 0.0

julia> Qn=Sn\Rn
3-element Vector{Float64}:
 -9.0
  0.5
  3.0

julia> ldiv!(Qn, Sn, Rn)
ERROR: MethodError: no method matching ldiv!(::Matrix{Int64}, ::Vector{Float64})

julia> SS=sparse([1 2 3;0 4 0;0 0 1])
3×3 SparseMatrixCSC{Int64, Int64} with 5 stored entries:
 1  2  3
 ⋅  4  ⋅
 ⋅  ⋅  1

julia> QS= spzeros(3)
3-element SparseVector{Float64, Int64} with 0 stored entries

julia> RS =sparse([1,2,3])
3-element SparseVector{Int64, Int64} with 3 stored entries:
  [1]  =  1
  [2]  =  2
  [3]  =  3

julia> QS=SS\RS
3-element Vector{Float64}:
 -9.0
  0.5
  3.0

julia> ldiv!(QS, SS, RS)
ERROR: MethodError: no method matching ldiv!(::SparseMatrixCSC{Int64, Int64}, ::SparseVector{Float64, Int64})
2

Have you read the docstring?

ldiv!(Y, A, B) -> Y

  Compute A \ B in-place and store the result in Y, returning the result.

  The argument A should not be a matrix. Rather, instead of matrices it should be a factorization object (e.g. produced by factorize or cholesky). The reason for this is that factorization itself is both
  expensive and typically allocates memory (although it can also be done in-place via, e.g., lu!), and performance-critical situations requiring ldiv! usually also require fine-grained control over the
  factorization of A.

  Examples
  ≡≡≡≡≡≡≡≡≡≡

  julia> A = [1 2.2 4; 3.1 0.2 3; 4 1 2];

  julia> X = [1; 2.5; 3];

  julia> Y = zero(X);

  julia> ldiv!(Y, qr(A), X);
4

The dosctring says:

ldiv!(Y, A, B)
(...)
The argument A should not be a matrix. Rather, instead of matrices it should be a factorization object (...)

Your code is doing:

julia> Sn=[1 2 3;0 4 0;0 0 1]
3×3 Matrix{Int64}:

and

julia> SS=sparse([1 2 3;0 4 0;0 0 1])
3×3 SparseMatrixCSC{Int64, Int64} with 5 stored entries:

as you can see from the type displayed after construction of your Sn/SS inputs, both are matrices, so directly contradict what’s in the docstring.

Keeping with your example:

julia> ldiv!(Qn, qr(Sn), Rn)
3-element Vector{Float64}:
 -9.0
  0.5
  3.0
1 Like
6

A rule of thumb is to only use sparse algorithms for matrices with less than 10% nonzero éléments. You have 45%

3 Likes
7

Thank you for your note. Aha, is this rule of thumb cited in the literature? Does this because of the complex structure of Sparse Matrix/vector than the normal matrix/vector (thus slower in performance if nonzero éléments is more than 10%)?

8

I quickly searched for the 10% rule but couldn’t find a reference to it. In any case, the benefits of sparse matrix code depends also of the operations you want to conduct. To make a sound decision, you should benchmark your code with dense or sparse matrices.
You are right, sparse matrix code is more complex and therefore takes more time if you have a lot on non-zero element in the matrix

1 Like
9

Thank you! I’ll keep that in my mind.

10

honestly 10% is very low sparsity. I would consider 1% closer to the breakeven point. Of course it depends on how big your matrices are and memory bandwidth vs compute bandwidth etc.

11

So, you recommend:
1- Sparsity of 50%∓1%?
2- With bigger matrices, the previous sparsity could be increased?