Loïc Grenié on Wed, 05 Feb 2014 11:16:51 +0100


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

Re: handling several substacks ?


2014-02-05 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.

    I think ellfacteur does something similar to what you are looking for.
  Look for alloc_scratch and the usage of the various tables allocated
  by alloc_scratch. I don't know whether there is some gerepile in that
  space. However, since your a[i] does not seem to use itself, you
  probably do not really need to gerepile.

       Hope this helps,

          Loïc