Pascal Molin on Wed, 05 Feb 2014 10:52:25 +0100

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

handling several substacks ?

I have to deal with the following kind of loops:

gploop(N=8000,K=400,L=8) = {
  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 */
      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" && /usr/local/bin/gcc-4.7 -o -bundle -flat_namespace -undefined suppress -O3 -Wall -fno-strict-aliasing -fomit-frame-pointer -fPIC stack.o "; -*-*/
#include <pari/pari.h>
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;

Description: Binary data