Nemocas / AbstractAlgebra.jl

Generic abstract algebra functionality in pure Julia (no C dependencies)
https://nemocas.github.io/AbstractAlgebra.jl/dev/index.html
Other
162 stars 62 forks source link

Memory allocation in diagonal_matrix #1108

Open a-kulkarn opened 2 years ago

a-kulkarn commented 2 years ago

Hello AbstractAlgebra folks,

Is this intended behaviour?

using Oscar
Qp = PadicField()
I = identity_matrix(QQ, 5)
I[1,1] === I[2,2] # evaluates to false

Ip = identity_matrix(Qp, 5)
Ip[1,1] === Ip[2,2] # evaluates to true

If so, I think it might be nice to make a note of this peculiarity in the documentation. In either case, I'm happy to follow this up with a PR.

tthsqe12 commented 2 years ago

~I think the second one is bad.~ Alright, the decision apparently was that both are acceptable.

thofma commented 2 years ago

No, this is intentional, see the discussion in https://github.com/Nemocas/AbstractAlgebra.jl/pull/692

thofma commented 2 years ago

It is documented as the first item here: http://nemocas.github.io/Nemo.jl/latest/developer/topics/#Aliasing-rules

tthsqe12 commented 2 years ago

hmm. I see this in there:

Constructors for polynomials, series and similar ring element objects (that are not matrices) that take an array as input, must ensure that the coefficients being placed into the object do not alias, even in part

Are we even doing this?

julia> using AbstractAlgebra

julia> R, x = PolynomialRing(ZZ, "x")
(Univariate Polynomial Ring in x over Integers, x)

julia> a = ZZ(1)
1

julia> p = R([a, a, a])
x^2 + x + 1

julia> coeff(p, 0) === coeff(p, 1)
true
thofma commented 2 years ago

Yes, you are constructing it (you are the "Constructor") and must ensure that the coefficients do not alias each other. We are not checking this :)

wbhart commented 2 years ago

I think it is badly worded. I think it means that users of such constructors...

tthsqe12 commented 2 years ago

Sure, so the documentation is meant to push the burden onto the user. But, of course the ownership problems don't end with this.

julia> R, x = PolynomialRing(ZZ, "x")
(Univariate Polynomial Ring in x over Integers, x)

julia> a = ZZ(1);

julia> array = [a, deepcopy(a), deepcopy(a)];
julia> p = R(array)
x^2 + x + 1

julia> array[2] += 1; p
x^2 + 2*x + 1
thofma commented 2 years ago

That is just the normal julia "no copy on construction" behavior which is unrelated to aliasing questions I believe.

I am happy to review/discuss any proposal to have a consistent and easy to use aliasing/mutation model for both matrices and polynomials or anything else. Bonus points if it can be implemented in julia.

a-kulkarn commented 2 years ago

Hmm...OK. The funny thing is that the "zero" entries in the identity matrix are unaliased, since it is created using zero_matrix.

The "Julian" solution might be to have a IdentityMatrix{T}, ZeroMatrix{T} <: MatElem type that implements the matrix interface, and for Nemo types to do their own thing with the identity_matrix, zero_matrix constructors, which it seems like they do anyways. This is yet-another-type-implementation, but I thought I would mention the idea.

Alternatively, one can introduce a matrix_malloc method. I'm not sure how often identity_matrix and zero_matrix are used as constructors, but it seems as though there is some active danger in doing so. Don't know if Tommy's comment from #692 is still relevant here.

The function works as long as you know that in your output matrix, none of the entries alias each other. At the moment the output of zero_matrix and identity_matrix satisfy this. We usually use those as the first argument to !-functions. But of course this works because of an implementation detail.

I am fully behind Nemocas/Nemo.jl#278, but maybe we need some more unsafe functions for matrices? Like add!! or something like that, where it is assumed that the entries of the output do not alias each other. I am not sure, but from our experience in the last 5 years I can tell that working with matrices is quite inefficient if we don't allow those entry-mutating functions.