handling several substacks ?

[Pascal Molin <molin@math.univ-paris-diderot.fr>]

I have to deal with the following kind of loops:

gploop(N=8000,K=400,L=8) = {

  my(m,a);

  a = vector(L,l,vector(K,k,exp(I*random(1.))));

  m = vector(L,l,vector(K,k,exp(I*random(1.))));

  for(n=1,N, /* complex outer loop, stack clean */

      my(i,j);

      i=random(L)+1;j=random(L)+1;

      for(k=1,K, /* trivial inner loop, stack demanding */

        a[i][k] = a[j][k] * m[i][k];

        );

      /* use A = a[i] here */

    );

}

gp handles them correctly, but I do not know how to translate this in libpari

using decent stack management. A gerepilecopy call on the vector a when

the stack gets empty is not affordable (this is impressive, see the c file attached, unless I have made a mistake).

And gaffect also makes copies.

I think I would like to allocate L+1 stack slices and make sure that each component

of some a[i] points only inside one of them, manually changing the stack pointer before

each k loop to a free slice.

This is the way mpc would work (thanks to Andreas last post, this gives a first solution).

Can this work and is this the right thing to do ? in this case I would be happy to have some hints on the way to do it (in particular the sizes).

-- 

Pascal Molin

/*-*- compile-command: "/usr/local/bin/gcc-4.7 -c -o stack.o -O3 -Wall -fno-strict-aliasing -fomit-frame-pointer -fPIC -I"/Users/pascal/devel/install/include" stack.gp.c && /usr/local/bin/gcc-4.7 -o stack.gp.so -bundle -flat_namespace -undefined suppress -O3 -Wall -fno-strict-aliasing -fomit-frame-pointer -fPIC stack.o "; -*-*/
#include <pari/pari.h>
/*
GP;install("cloop","D50000,L,D1000,L,D8,L,p","cloop","./stack.so");
*/
GEN
cloop(long N, long K, long L, long prec)
{
  pari_sp ltop = avma;
  GEN a, m;
  {
    /* init */
    long l;
    a = cgetg(L+1, t_VEC);
    m = cgetg(L+1, t_VEC);
    for (l = 1; l <= L; ++l)
    {
      long k;
      GEN p1, p2;
      p1 = cgetg(K+1, t_VEC);
      p2 = cgetg(K+1, t_VEC);
      for (k = 1; k <= K; ++k)
      {
          gel(p1, k) = gexp(gmul(gen_I(), genrand(real_1(prec))), prec);
          gel(p2, k) = gexp(gmul(gen_I(), genrand(real_1(prec))), prec);
      }
      gel(a, l) = p1;
      gel(m, l) = p2;
    }
  }
  {
    /* main loop */
    pari_sp av1 = avma, lim1 = stack_lim(avma, 1);
    long n;
    for (n = 1; n <= N; ++n)
    {
        long i, j, k;
        i = random_Fl(L)+1;
        j = random_Fl(L)+1;
          for (k = 1; k <= K; ++k)
          {
            gel(gel(a, i), k) = gmul( gel(gel(a, j), k), gel(gel(m, i), k));
            if (avma > lim1)
              a = gerepilecopy(av1, a);
      }
    }
  }
  /* use A = a[i] here */
  avma = ltop;
  return gen_0;
}

Attachment: stack.gp
Description: Binary data