davidavdav / InplaceLinalg.jl

Julia macro package supporting in-place linear algebra
MIT License
0 stars 0 forks source link

@inplace should not scale aguments #3

Open bsxfan opened 5 years ago

bsxfan commented 5 years ago

Consider:

julia> @macroexpand @inplace C -= R*2*S
:(InplaceLinalg.C_AB!(C, 1, -R, 2, S))

If R is a matrix (or vector) then -R is not done inplace---and unnecesarily in any case, since gemm! will be doing the scaling anyway. The macro should instead do:

julia> @macroexpand @inplace C -= R*2*S
:(InplaceLinalg.add_update!(C, 1, -,R, 2, S))
bsxfan commented 5 years ago

The problem also happpens with multiplicative updates, e.g.:

julia> @macroexpand @inplace B = -A*B
:(InplaceLinalg.mult_update!(B, -A, Val(2)))

This should do instead:

mult_update!(B, -1, A, Val(3))

And this has the same problem:

julia> @macroexpand @inplace B = -2A*B
:(InplaceLinalg.mult_update!(B, -2A, Val(2)))

and this:

julia> @macroexpand @inplace B = 2A*B
:(InplaceLinalg.mult_update!(B, 2A, Val(2)))

A probably more difficult case is:

julia> @macroexpand @inplace B = -B*A
:(InplaceLinalg.C_AB!(B, 0, 1, -B, A))

because it fails to recognize that this can be a multiplicative update.

davidavdav commented 5 years ago

OK, I got the cases

@inplace B = -A*B
@inplace B = -B*A

which were relatively simple. But I've spent too much time now trying to write a function that simply flattens out multiplications.

davidavdav commented 5 years ago

Some even more esoteric cases aren't handled properly, like

@inplace B = -2A*-B
:(InplaceLinalg.mult_update!(B, -1, -2, A, Val(4)))

whereas things like

@macroexpand @inplace B = -A*-B
:(InplaceLinalg.mult_update!(B, A, Val(2)))

which is intentional.

I suppose I could make checksigns also be on the watchout for negative numbers. Currently it only finds - operators.

davidavdav commented 5 years ago

Yes, I can do this but that is shifting the problem somewhere else, and I can try to start to fix that, but really, I think it is not wise if I start to look for numbers. Because I can't start looking for scalars.

So the question then is if we want to parse B = -B*A and B = -2B*A, do we also need to be able to parseB = (-a*B)*A, which is the generalization from -2 to -a?

bsxfan commented 5 years ago

Commit 0c033f9, fixes the problem of scaling the sedond term in additive updates. The first example mentioned in this issue, C -= R*2*S is now working properly.

The same commit now also handles additive updates with differences of terms, e.g.:

C = 2C - A*B
C -= 3C - A*B