RFC - Change signature of methods A_ldiv_B!(A, B) to A_ldiv_B!(X, A, B)

[KristofferC]

Currently, the methods for solving linear system of equations in place (A_ldiv_B!) takes two arguments A and B. In some cases, calling A_ldiv_B!(A, B) will update B, in some cases not, see: #201.

There are many solvers that does not support updating B directly but instead needs a preallocated solution array passed into them to store the solution in. This is true for i.e. UMFPACK which means that currently there is no way to call A_ldiv_B! for sparse A and have it not allocate the solution array at each call.

It also means that there is no clean way to keep your RHS B intact while calling in place solvers. Instead you need to copy B before the call since B gets overwritten.

I believe that a more natural signature for these methods would be A_ldiv_B!(X, A, B) where X is the preallocated array to store the solution. This has symmetry with A_mul_B!(X, A, B) and it leaves the RHS untouched.

I might have missed a good reason why the current signatures are the way there are but from my point of view this looks cleaner. It is also the path I took in https://github.com/KristofferC/Pardiso.jl where I have the following solve methods:

• solve(ps, A, B) solves AX=B and returns X
• solve!(ps, X, A, B) solves AX=B and stores it in X

Any comments on the sensibility and feasibility for this suggestion?

[timholy]

I really love this. It's more convenient to use and more consistent.

[ViralBShah]

+1

[lindahua]

Sometimes, we may want to solve
XA = B. So there is distinction between left and right division.

[KristofferC]

@lindahua Are you referring to the solve functions I gave as an example or do you mean the Base functions? From what I can see, there are currently no in place right division in Base.

[simonbyrne]

+1

[simonbyrne]

@lindahua: I think @KristofferC is proposing adding additional A_ldiv_B! methods, not a generic solve! method.

[lindahua]

Thanks for clarification. +1

[tkelman]

The consistency here sounds great to me. We'd want to be careful about aliasing scenarios, I imagine A_ldiv_B!(B, A, B) might end up commonly used for types where that makes sense.

[simonbyrne]

I imagine A_ldiv_B!(B, A, B) might end up commonly used for types where that makes sense.

I don't think this is a pattern we should encourage, as it is going to break in some cases, e.g.:

julia> X = randn(5,5)
5x5 Array{Float64,2}:
  1.33502    0.168431  -1.24184     0.152598   0.368737
 -0.984189  -0.112196  -0.0336856  -0.281819   0.881885
  1.35078    0.7121    -1.2564      0.316167  -0.14738 
  1.01016   -1.84001   -0.577902    0.286557  -1.2575  
  0.560689   2.57815    1.20324    -1.54285    1.1421  

julia> A_mul_B!(X,X,X)
5x5 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

As a blatant plug, can I suggest the Inplace{N} type I use in InplaceOps.jl?

[KristofferC]

Regarding A_ldiv_B!(B, A, B) isn't that in the "consenting adults" territory? As @simonbyrne pointed out, the same issue already exists for in place matrix multiplication.

[toivoh]

It seems to me that the problem with this kind of aliasing situation is that depending on the type of A, it might be fine or it might not. That feels more like a trap.

[KristofferC]

To be honest, I don't really see the issue with the potential aliasing problem. I would argue if you are using A_ldiv_B! it means you are an advanced user and should know when it is safe or not. We do exactly the same now with A_mul_Bx!.

[andreasnoack]

I think the better approach is that of BLAS and LAPACK where you only use output arrays if the algorithm requires. Therefore, it would be preferable if Pardiso and SuiteSparse changed their signatures to update b in place. However, that might be too much to ask for so I'd be okay with adding methods with dedicated output arrays for generic programming purposes, but we should keep the in-place versions.

I completely agree with @KristofferC's last comment that users of the ! should be expected to know the algorithms or read the documentation to avoid aliasing problems.

[tkelman]

users of the ! should be expected to know the algorithms or read the documentation to avoid aliasing problems.

Writers of the ! methods will have to write that documentation to spell out when it is or isn't safe for the second option to work :)

[andreasnoack]

I think the approach should be to note that aliasing of input arguments is illegal in general and comment specifically when aliasing is ok.

[simonbyrne]

I would suggest that the user should operate on the assumption that the 3-argument form can't be aliased, and we should provide the 2-argument form (or use an Inplace{2} marker) for when it is okay: then if the method can't operate inplace, the user sees a MethodError, rather than an obtaining incorrect result.

[andreasnoack]

The easy thing about the BLAS/LAPACK approach is that aliasing is never okay.

I'd also like to add that Julia/LLVM has a hard time optimizing the loads and stores of arrays because of the aliasing. This only applies to the functions implemented in Julia but anyway. The more vector arguments, the more potential aliasing to make optimizations harder.

[tkelman]

Bump. It seemed like almost everyone was in favor of this.

[KristofferC]

I've started on UMFPACK, at least for real types.

[andreasnoack]

I just looked at CHOLMOD's cholmod_solve2 and there you'll have to supply a handle to the result vector, i.e. a Ptr{Ptr{T}}. It that the same for UMFPACK? We'd need to make sure that the SuiteSparse is not messing around with Julia's allocations.

[JaredCrean2]

This is great. Looking forward to when this is merged.

[KristofferC]

UMFPACK is much easier. Just swap a similar to an inplace function and add a solve(lu::UmfpackLU{Float64,$itype}, b::VecOrMat{Float64}, typ::Integer) = solve!(similar(b), lu, b)

[KristofferC]

For the LAPACK versions that solve AX = B and stores the solution in B do we just write methods similar to described at: JuliaLang/julia#16433 (comment)

[andreasnoack]

Are you asking if we should support three argument A_ldiv_B! methods for the LAPACK types?

[KristofferC]

Yes.

[andreasnoack]

For now, we should probably support the three argument version for all A_ldiv_B! such that we have the same API for dense and sparse matrices. Personally, I'd prefer that SuiteSparse fixed the API to match LAPACK.

[stevengj]

3-arg form added in JuliaLang/julia#16702

[StefanKarpinski]

Can this be closed then?

[KristofferC]

Still lacks CHOLMOD if I understand things correctly.

[mfalt]

I made an implementation for A_ldiv_B! using CHOLMOD solve2 here: https://github.com/mfalt/CholmodSolve2.jl

It currently only supports Float64 (not Complex128), but that should be (relatively) easy to add. I solved the problem of references by making sure CHOLMOD does all the allocations. I think the code works as expected, I have not found any memory leaks at least.
I am not sure how it will work in parallel settings, and the work space will not be free'd until julia quits. Maybe you have some input on this?

There were some problems that makes the code a bit ugly (and unsafe?).
I think I found a "bug" in solve2 which I emailed the creator about. The problem is that although workspace variables are submitted of the right size, CHOLMOD will reallocate them if they are not of right dimensions. This could be all right if CHOLMOD didn't mess up the dimensions on return, making it impossible to not allocate on subsequent calls using the standard interface.

For example if RHS has less than 4 columns (say 1), CHOLMOD still wants to allocate work-space Y to be 4-by-n. If less is supplied (or more) it will free and reallocate 4-by-n workspace.
If Y is supplied as 4-by-n (i.e Y->nzmax = 4*n), then the check in CHOLMOD that Y is 4-by-n passes fine.
But a couple of lines later, when transposing the output, CHOLMOD wants Y to be 1-by-n.
This is simply solved by setting Y->nrow = 1; Y->d = 1 in CHOLMOD.
However, Y is returned in this state(!), so the next call to solve2 doesn't recognize that Y->nzmax = 4*n and will free and then reallocate Y of the exact same size.

I solved this by manually setting Y->ncols, Y->nrows, Y->d from julia, if Y->nzmax is large enough. So the current implementation will do O(1) allocations, as long as the current problem is not larger than any previous problem.

[ViralBShah]

I believe this can be closed but please do reopen if necessary.