question on matrices sharing some elements

[Giacomo Mulas <gmulas@oa-cagliari.inaf.it>]

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 need to check for the existence of solutions of a system of diophantine 
equations, which I need to repeat a large number of times with the
coefficients matrix constant and only the constant terms changing.
At the top I set

  pari_sp ltop = avma;

Then I begin by creating the coefficients matrix with

GEN coeffs1 = zeromatcopy(lines, cols);

Then I set up the nonzero coefficients with a bunch of

  gmael2(coeffs1,i,j) = stoi(1);

and compute the SNF of coeffs1 with

  GEN snf1 = ZM_snf(coeffs1);

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);
  }
  gel(coeffs2,cols+1) = zerocol(lines);

The intended meaning is that the extended matrix is identical to the
coefficients matrix, just with the added column of constant terms, so
identical columns are not duplicated, but pointed to by both matrices.
The additional column is a new column vector.

Then I set up the nonzero constant terms with a bunch of

  gmael2(coeffs2,i,cols+1) = stoi(1);

At this stage, however, somethng is wrong with the coefficients matrix, 
because if I add a

  pari_printf("coeffs1 after gc:\n%Ps\n", coeffs1);

some of its terms were modified.

Would some PARI/GP guru help me to understand what I did wrong here?

A skeletal C source code follows to demonstrate the problem.

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? Or should I
explicitely enumerate all the leaves of a tree I want to preserve, when
doing garbage collection?

Thanks in advance for any help, bye
Giacomo


#include <pari/pari.h>

void main() {
  long lines = 4;
  long cols1 = 3;
  long cols2 = 4;
  long i,j;

  pari_init(100000000, 0);

  pari_sp ltop = avma;

  /* define the matrix of coefficients */
  GEN coeffs1 = zeromatcopy(lines, cols1);

  gmael2(coeffs1,2,1) = stoi(1);
  gmael2(coeffs1,3,1) = stoi(-1);
  gmael2(coeffs1,1,2) = stoi(1);
  gmael2(coeffs1,2,2) = stoi(-1);
  gmael2(coeffs1,1,3) = stoi(1);
  gmael2(coeffs1,3,3) = stoi(1);

  /* compute the SNF */

  GEN snf1 = ZM_snf(coeffs1);

  pari_printf("coeffs1:\n%Ps\n", coeffs1);
  pari_printf("snf1:\n%Ps\n\n", snf1);

  /* define the extended matrix */

  GEN coeffs2 = cgetg(cols2+1, t_MAT);
  for (j=1; j<=cols1; j++) {
    gel(coeffs2,j) = gel(coeffs1,j);
  }
  gel(coeffs2,cols2) = zerocol(lines);

  gmael2(coeffs2,2,4) = stoi(2);

  /* now only the extended matrix should be affected, instead... */
  pari_printf("coeffs1 changed!:\n%Ps\n", coeffs1);
  /* but the extended matrix is as expected */
  pari_printf("coeffs2 :\n%Ps\n", coeffs1);

  /* compute the SNF of the extended matrix */
  GEN snf2 = ZM_snf(coeffs2);

  pari_printf("snf2:\n%Ps\n", snf1);
  /* nonzero terms of snf1 and snf2 are identical, this system has solutions
  */

  /* I don't need coeefs1 anymore, do GC */
  gerepileall(ltop, 3, &coeffs2, &snf1, &snf2);

  /* change the constant terms, compute another SNF of the extended matrix*/
  gmael2(coeffs2,3,4) = stoi(1);
  GEN snf3 = ZM_snf(coeffs2);

  pari_printf("snf3 :\n%Ps\n", snf3);
  /* this time nonzero terms of snf1 and snf3 are not identical, this system
   * has no solutions */

  /* do GC */
  gerepileall(ltop, 4, &coeffs2, &snf1, &snf2, &snf3);

  pari_close();

  return ;

}


--
_________________________________________________________________

Giacomo Mulas <gmulas@oa-cagliari.inaf.it>
_________________________________________________________________

INAF - Osservatorio Astronomico di Cagliari
via della scienza 5 - 09047 Selargius (CA)

tel.   +39 070 71180244
mob. : +39 329  6603810
_________________________________________________________________

"When the storms are raging around you, stay right where you are"
                         (Freddy Mercury)
_________________________________________________________________