Bill Allombert on Thu, 29 Sep 2005 16:38:10 +0200


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

Re: PARI and POSIX threads


On Thu, Aug 18, 2005 at 05:04:09PM +0200, Bill Allombert wrote:
> On Mon, Aug 15, 2005 at 11:24:58PM +0200, Philippe Elbaz-Vincent wrote:
> > Hi,
> > 
> > Is there a solution to this problem ?
> > [is it technically doable to run several instances of PARI inside the same 
> > program ?]
> 
> I am not sure the answer is known.
> 
> What you can try is to mark avma, bot and top 'thread-local' using the
> '__thread' keyword (assuming this is supported) and use
> allocatemoremem() or switch_stack() to allocate a separate stack in each
> threads.  To copy objects between threads, clone them.
> 
> We would be very interested to know if that can be made to work.

Since the best way is to try, here a patch that seems to allow to run
concurrent PARI threads assuming thread-local storage is supported:

This:
1) add a THREAD macro that default to the empty string
2) mark avma,top and bot with THREAD (as a storage class).
3) add 2 public functions pari_thread_init() and pari_thread_close().

As such the patch does absolutly nothing until you compile PARI
with -DTHREAD=__thread on a system that support thread-local storage.
(You probably also want -D_REENTRANT).
The simplest is to do:
make clean
make gp CC_FLAVOR=-DTHREAD=__thread
make install

Once you have installed the resulting library, you can try
the attached example program thread.c.

Here a commented version:

#include <pari/pari.h>
#include <pthread.h>

void * mydet(void * arg)
{
  GEN M=(GEN)arg;
  GEN F;
  pari_thread_init(8000000); // Initialize a local PARI stack for this
                             // thread.
  F=gclone(det(M));          // compute the determinant using the local
                             // stack
                             // clone the result so we can free the
                             // stack
  pari_thread_close();       // Free the PARI stack
  pthread_exit((void *)F);   // End the thread and return F.
  return F;                  // spurious
}

void * myfactor(void * arg)  //exactly the same as above.
{
  GEN N=(GEN)arg;
  GEN F;
  pari_thread_init(4000000);
  F=gclone(factor(N));
  pari_thread_close();
  pthread_exit((void *)F);
  return F;
}

int main(void)
{
  GEN M,N1,N2, F1,F2,D;
  pthread_t th1, th2, th3;
  pari_init(4000000,500000);         //Initialise the main PARI stack
                                     //and global objects (gen_0, etc.)
                                     //We rely on the assumption that 
                                     //the threads will not modify any
                                     //non-constant global objects.
                                     //(Pi,log2,etc.)
  N1=addis(gpowgs(gen_2,256),1);     //do computation in the main PARI
                                     //stack
  N2=subis(gpowgs(gen_2,193),1);
  M=mathilbert(80);
  pthread_create(&th1,NULL,&myfactor,(void*)N1);//Start threads
  pthread_create(&th2,NULL,&myfactor,(void*)N2);
  pthread_create(&th3,NULL,&mydet,(void*)M);
  pthread_join(th1,(void*)&F1);//Wait until they terminate and get
                               //theirs results
  pthread_join(th2,(void*)&F2);
  pthread_join(th3,(void*)&D);
  pariputsf("F1=%Z\nF2=%Z\nlog(D)=%Z\n",F1,F2,glog(D,3));
                                            //display the results.
  return 0;
}

Example run on a quadri-opteron:
$ gcc -I amd64/include -L amd64/lib -lpari -lgmp thread.c -O3 -othread -Wall -lpthread
F1=[1238926361552897, 1; 93461639715357977769163558199606896584051237541638188580280321, 1]
F2=[13821503, 1; 61654440233248340616559, 1; 14732265321145317331353282383, 1]
log(D)=-8726.787756120586376
LD_LIBRARY_PATH=amd64/lib ./thread  12,18s user 0,16s system 264% cpu 4,668 total

So we used 4,668s instead of 12,18s with a single thread, and more
importantly we did not randomly crash and even the output is correct.

Cheers,
Bill
? amd64
? patch
? thread
? thread.c
Index: src/headers/paridecl.h
===================================================================
RCS file: /home/cvs/pari/src/headers/paridecl.h,v
retrieving revision 1.513
diff -u -r1.513 paridecl.h
--- src/headers/paridecl.h	26 Sep 2005 15:00:15 -0000	1.513
+++ src/headers/paridecl.h	29 Sep 2005 13:42:42 -0000
@@ -1198,6 +1198,8 @@
 void    err_leave(void **v);
 GEN     forcecopy(GEN x);
 void    freeall(void);
+void    pari_thread_init(size_t parisize);
+void    pari_thread_close(void);
 GEN     gcopy(GEN x);
 GEN     gcopy_i(GEN x, long lx);
 GEN     gerepile(pari_sp ltop, pari_sp lbot, GEN q);
Index: src/headers/paristio.h
===================================================================
RCS file: /home/cvs/pari/src/headers/paristio.h,v
retrieving revision 1.30
diff -u -r1.30 paristio.h
--- src/headers/paristio.h	6 Jul 2005 19:12:32 -0000	1.30
+++ src/headers/paristio.h	29 Sep 2005 13:42:42 -0000
@@ -74,7 +74,7 @@
 #define TEXSTYLE_PAREN	2
 #define TEXSTYLE_BREAK	4
 
-extern pari_sp avma,bot,top;
+extern pari_sp THREAD avma,bot,top;
 #define DISABLE_MEMUSED (size_t)-1
 extern size_t memused;
 extern byteptr diffptr;
Index: src/headers/parisys.h
===================================================================
RCS file: /home/cvs/pari/src/headers/parisys.h,v
retrieving revision 1.10
diff -u -r1.10 parisys.h
--- src/headers/parisys.h	1 Sep 2005 17:02:43 -0000	1.10
+++ src/headers/parisys.h	29 Sep 2005 13:42:42 -0000
@@ -57,6 +57,9 @@
 #  define INLINE_IS_STATIC
 #  define INLINE static
 #endif
+#ifndef THREAD
+#  define THREAD
+#endif
 
 #if defined(_WIN32) || defined(__CYGWIN32__)
 /* ANSI C does not allow to longjmp() out of a signal handler, in particular,
Index: src/language/init.c
===================================================================
RCS file: /home/cvs/pari/src/language/init.c,v
retrieving revision 1.272
diff -u -r1.272 init.c
--- src/language/init.c	3 Sep 2005 16:36:12 -0000	1.272
+++ src/language/init.c	29 Sep 2005 13:42:42 -0000
@@ -48,7 +48,7 @@
 ulong   DEBUGFILES, DEBUGLEVEL, DEBUGMEM, compatible;
 ulong   prec, precdl;
 ulong   init_opts = INIT_JMPm | INIT_SIGm;
-pari_sp bot = 0, top = 0, avma;
+THREAD  pari_sp bot = 0, top = 0, avma;
 size_t memused;
 
 gp_data *GP_DATA = NULL;
@@ -638,6 +638,18 @@
 {
   return gp_init_entrees(new_fun_set? pari_modules: pari_oldmodules,
                          functions_hash, force);
+}
+
+void
+pari_thread_init(size_t parisize)
+{
+  init_stack(parisize);
+}
+
+void
+pari_thread_close(void)
+{
+  free((void *)bot);
 }
 
 /* initialize PARI data. You can call pari_addfunctions() first to add other
#include <pari/pari.h>
#include <pthread.h>

void * mydet(void * arg)
{
  GEN M=(GEN)arg;
  GEN F;
  pari_thread_init(8000000);
  F=gclone(det(M));
  pari_thread_close();
  pthread_exit((void *)F);
  return F;
}

void * myfactor(void * arg)
{
  GEN N=(GEN)arg;
  GEN F;
  pari_thread_init(4000000);
  F=gclone(factor(N));
  pari_thread_close();
  pthread_exit((void *)F);
  return F;
}

int main(void)
{
  GEN M,N1,N2, F1,F2,D;
  pthread_t th1, th2, th3;
  pari_init(4000000,500000);
  N1=addis(gpowgs(gen_2,256),1);
  N2=subis(gpowgs(gen_2,193),1);
  M=mathilbert(80);
  pthread_create(&th1,NULL,&myfactor,(void*)N1);
  pthread_create(&th2,NULL,&myfactor,(void*)N2);
  pthread_create(&th3,NULL,&mydet,(void*)M);
  pthread_join(th1,(void*)&F1);
  pthread_join(th2,(void*)&F2);
  pthread_join(th3,(void*)&D);
  pariputsf("F1=%Z\nF2=%Z\nlog(D)=%Z\n",F1,F2,glog(D,3));
  return 0;
}