Giacomo Mulas on Wed, 26 Mar 2014 13:54:35 +0100 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
question on matrices sharing some elements |
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) _________________________________________________________________