Bill Allombert on Wed, 26 Mar 2014 16:12:41 +0100

[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

Re: question on matrices sharing some elements

On Wed, Mar 26, 2014 at 01:54:05PM +0100, Giacomo Mulas wrote:
> Hello, I am a new (ab)user of pari, please CC to me any replies, as I am not
> (yet) subscribed to the list.
> I only want to use the library, and I need only a very small part of its
> functionality, namely computing the Smith normal form of an integer matrix.
> First off, therefore, if anyone can point me to a simpler library doing just
> this, in particular if it is optimised for sparse matrices containing only
> small numbers (the sum of the absolute values of any column is >=4), this
> would be most welcome.

I would suggest linbox <> which should be faster but I
doubt this is simpler to use.

> Then I create the extended matrix with
>   GEN coeffs2 = cgetg(cols+2, t_MAT);
>   for (j=1; j<=cols; j++) {
>     gel(coeffs2,j) = gel(coeffs1,j);
>   }

This is where the bug lies: you need to copy gel(coeffs1,j):
gel(coeffs2,j) = gcopy(gel(coeffs1,j));
or even simply
GEN coeffs2 = gcopy(coeffs1);
otherwise the columns of coeffs2 and coeffs1 are sharing the same memory.

>   /* but the extended matrix is as expected */
>   pari_printf("coeffs2 :\n%Ps\n", coeffs1);

you meant to do
pari_printf("coeffs2 :\n%Ps\n", coeffs2);

> I also have another question: my understanding of gerepileall is that it
> preserves recursively also all the leaves of a tree object, is this correct?
> If I write:
> gerepileall(ltop, 1, &coeffs2);
> and coeffs2 is a t_MAT object, the stack will be cleaned but preserving
> coeffs2, its columns and all the elements of each column, right?

Yes, this works recursively.