Changes from Dylan Stark to allow registration of topology

routines and introducing a physical to logical processor
mapping.

Now it is possible to register custom topology routines if
you don't like the two ones we have at the moment (the
traditional one and the new one from Frank Loeffler).

You can also provide routines to optimise the placement of
neighbouring processes on the processors by giving a custom
mapping from the numbers PUGH uses for the processors to
the actual MPI processor.

Please treat this registration as semi-volatile for the moment
and coordinate with Dylan if you fancy using it.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusPUGH/PUGH/trunk@475 b61c5cb5-eaca-4651-9a7a-d64986f99364

10 files changed, 1058 insertions, 303 deletions

diff --git a/interface.ccl b/interface.ccl
index a445a0f..95c3b2a 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -2,3 +2,7 @@
 # $Header$
 
 implements: Driver
+
+INCLUDE HEADER: Register.h in pugh_Register.h
+
+USES INCLUDE: pugh_Register.h
diff --git a/param.ccl b/param.ccl
index 5a751ca..da5af5b 100644
--- a/param.ccl
+++ b/param.ccl
@@ -112,6 +112,14 @@ BOOLEAN enable_all_storage "Enable storage for all GFs?"
 
 # Processor topology
 
+restricted:
+
+KEYWORD physical2logical "Physical process to logical process mapping method to use" STEERABLE = RECOVER
+{
+  "direct"  :: "Maps MPI IDs directly to IJKs"
+  "example" :: "Maps MPI IDs directly to IJKs using a lookup table"
+} "direct"
+
 KEYWORD processor_topology "How to determine the processor topology" STEERABLE = RECOVER
 {
   "manual"        :: "Specified by proc_top_nx etc"
@@ -119,6 +127,8 @@ KEYWORD processor_topology "How to determine the processor topology" STEERABLE =
   "automatic_old" :: "Automatically generated (old method)"
 } "automatic"
 
+private:
+
 INT processor_topology_1d_x "No of Procs in X direction" STEERABLE = RECOVER
 {
   0:*  :: "See proc_topology"
diff --git a/schedule.ccl b/schedule.ccl
index 3b4e2c2..26a11da 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -6,6 +6,16 @@ schedule PUGH_Startup at STARTUP as Driver_Startup
   LANG:C
 } "Startup routine"
 
+schedule PUGH_RegisterPUGHP2LRoutines at STARTUP
+{
+  LANG:C
+} "Register Physical to Logical process mapping routines"
+
+schedule PUGH_RegisterPUGHTopologyRoutines at STARTUP
+{
+  LANG:C
+} "Register topology generation routines routines"
+
 schedule PUGH_Report at BASEGRID
 {
   LANG:C
diff --git a/src/P2LMappings.c b/src/P2LMappings.c
new file mode 100644
index 0000000..85461fd
--- /dev/null
+++ b/src/P2LMappings.c
@@ -0,0 +1,269 @@
+ /*@@
+   @file      P2LMappings.c
+   @date      Sun Oct 30 17:02:58 CST 2005
+   @author    Dylan Stark
+   @desc
+              Provides routines for mapping between physical and
+              logical numberings of processors.
+   @enddesc
+   @version $Header$
+ @@*/
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "cctk.h"
+
+#include "pugh_Register.h"
+
+static const char *rcsid = "$Header$";
+
+CCTK_FILEVERSION(CactusPUGH_PUGH_P2L_c);
+
+
+/********************************************************************
+ *********************     Local Data Types   ***********************
+ ********************************************************************/
+
+/********************************************************************
+ ********************* Local Routine Prototypes *********************
+ ********************************************************************/
+
+static void *CreateDirectTranslation(int dim, int *nprocs, int np);
+static int   DestroyDirectTranslation(void *table);
+static int   DirectPhysicalToLogical(void *table, int proc);
+static int   DirectLogicalToPhysical(void *table, int ijk);
+
+static void *CreateExampleTranslation(int dim, int *nprocs, int np);
+static int   DestroyExampleTranslation(void *table);
+static int   ExamplePhysicalToLogical(void *table, int proc);
+static int   ExampleLogicalToPhysical(void *table, int ijk);
+
+/********************************************************************
+ ********************* Other Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ *********************     Local Data   *****************************
+ ********************************************************************/
+
+/********************************************************************
+ *********************  Scheduled Routines  *************************
+ ********************************************************************/
+/* Fortran wrappers appear immediately after C definition of function */
+
+ /*@@
+   @routine    PUGH_RegisterDefaultRoutines
+   @date       Wed Dec 14 16:51:22 2005
+   @author     Tom Goodale
+   @desc 
+   Registers the default translation and topology generation routines.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+ 
+ @@*/
+
+int PUGH_RegisterPUGHP2LRoutines(void)
+{
+  PUGH_RegisterP2L(CreateDirectTranslation,
+                   DestroyDirectTranslation,
+                   DirectPhysicalToLogical,
+                   DirectLogicalToPhysical,
+                   "direct");
+
+  PUGH_RegisterP2L(CreateExampleTranslation,
+                   DestroyExampleTranslation,
+                   ExamplePhysicalToLogical,
+                   ExampleLogicalToPhysical,
+                   "example");
+
+  return 0;
+}
+
+/********************************************************************
+ *********************     External Routines   **********************
+ ********************************************************************/
+
+
+/********************************************************************
+ *********************  Local Routines  *****************************
+ ********************************************************************/
+
+ /*@@
+   @routine    CreateDirectTranslation
+   @date       Wed Dec 14 19:00:57 2005
+   @author     Tom Goodale
+   @desc 
+   Creates a dummy translation table as one isn't necessary for
+   the direct case.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+ 
+ @@*/
+void *CreateDirectTranslation(int dim, int *nprocs, int np)
+{
+  static int dummy = 1;
+  return (void *)&dummy;
+}
+
+ /*@@
+   @routine    DestroyDirectTranslation
+   @date       Wed Dec 14 19:00:57 2005
+   @author     Tom Goodale
+   @desc 
+   Routine which would destroy the translation table if there was one.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+ 
+ @@*/
+static int DestroyDirectTranslation(void *table)
+{
+  return 0;
+}
+
+ /*@@
+   @routine    DirectPhysicalToLogical
+   @date       Wed Dec 14 19:02:56 2005
+   @author     Tom Goodale
+   @desc 
+   Direct map of the physical processor number to the logical one.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+ 
+ @@*/
+static int DirectPhysicalToLogical(void *table, int proc)
+{
+  return proc;
+}
+
+ /*@@
+   @routine    DirectLogicalToPhysical
+   @date       Wed Dec 14 19:02:56 2005
+   @author     Tom Goodale
+   @desc 
+   Direct map of the logical processor number to the physical one.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+ 
+ @@*/
+static int DirectLogicalToPhysical(void *table, int ijk)
+{
+  return ijk;
+}
+
+/*@@
+   @routine    CreateExampleTranslation
+   @date       Sat Oct  8 22:24:53 CDT 2005
+   @author     Dylan Stark
+   @desc
+               Allocates an int[2][np] array where the first
+               row maps id->ijk and the second row maps ijk->id.
+   @enddesc
+ @@*/
+void *CreateExampleTranslation(int dim, int *nprocs, int np)
+{
+  int i;
+  int **translation;
+  
+  translation    = (int **)malloc(2*sizeof(int *));
+  if (translation)
+  {
+    translation[0] = malloc(np*sizeof(int));
+    translation[1] = malloc(np*sizeof(int));
+
+    if (translation[0] && translation[1])
+    {
+      for (i=0; i<np; i++)
+      {
+        translation[0][i]=i;
+        translation[1][i]=i;
+      }
+    }
+    else
+    {
+      if (translation[0])
+      {
+        free(translation[0]);
+      }
+      if (translation[1])
+      {
+        free(translation[1]);
+      }
+      
+      free(translation);
+      translation = NULL;
+    }
+  }
+
+  return (void *)translation;
+}
+
+/*@@
+   @routine    DestroyExampleTranslation
+   @date       Sat Oct  8 22:31:52 CDT 2005
+   @author     Dylan Stark
+   @desc
+               Deallocate memory for the translation.
+   @enddesc
+ @@*/
+static int DestroyExampleTranslation(void *table)
+{
+  int **translation = (int **)table;
+
+  free(translation[0]);
+  free(translation[1]);
+  free(translation);
+
+  return 0;
+}
+
+
+/*@@
+   @routine    ExamplePhysicalToLogical
+   @date       Sat Oct  8 21:45:52 CDT 2005
+   @author     Dylan Stark
+   @desc
+     Take an MPI processor number and map it to some logical index.
+   @enddesc
+ @@*/
+static int ExamplePhysicalToLogical(void *table, int proc)
+{
+  int **translation = (int **)table;
+  
+  return translation[0][proc];
+}
+
+ /*@@
+   @routine    LogicalToPhysical
+   @date       Sat Oct  8 21:48:52 CDT 2005
+   @author     Dylan Stark
+   @desc
+     Take a logical index and map it to a MPI processor number.
+   @enddesc
+@@*/
+static int ExampleLogicalToPhysical(void *table, int ijk)
+{
+  int **translation = (int **)table;
+  
+  return translation[1][ijk];
+}
diff --git a/src/Register.h b/src/Register.h
new file mode 100644
index 0000000..5240375
--- /dev/null
+++ b/src/Register.h
@@ -0,0 +1,42 @@
+ /*@@
+   @header    Register.h
+   @date      Wed Dec 14 17:43:42 2005
+   @author    Tom Goodale
+   @desc 
+   This is a template header function
+   @enddesc
+   @version $Header$
+ @@*/
+
+#ifndef _PUGH_REGISTER_H_
+#define _PUGH_REGISTER_H_ 1
+
+#ifdef __cplusplus
+extern "C" 
+{
+#endif
+
+int PUGH_RegisterGenerateTopology(int (*GenerateTopology) (int dim,
+                                                           int total_procs, 
+                                                           const int *nsize, 
+                                                           const int *nghostzones, 
+                                                           int *nprocs), 
+                                  const char *topologyname);
+
+int PUGH_RegisterP2L(
+                     void *(*CreateP2LTranslation) (int dim, 
+                                                    int *nprocs, 
+                                                    int np),
+                     int   (*DestroyP2LTranslation) (void *table),
+                     int   (*PhysicalToLogical) (void *table, 
+                                                 int proc),
+                     int   (*LogicalToPhysical) (void *table, 
+                                                 int num),
+                     const char *p2lname);
+
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* _PUGH_REGISTER_H_ */
diff --git a/src/Registration.c b/src/Registration.c
new file mode 100644
index 0000000..549eff7
--- /dev/null
+++ b/src/Registration.c
@@ -0,0 +1,312 @@
+/*@@
+  @file      Registration.c
+  @date      Sat Oct  8 19:29:25 CDT 2005
+  @author    Dylan Stark (TommyD)
+  @desc
+             Allows registration of topology and 
+             physical to logical processor mapping routines.
+
+             Mostly taken from work done by Thomas Dramlitsch.
+  @enddesc
+  @version   $Header$
+ @@*/
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "Register.h"
+
+static const char *rcsid = "$Header$";
+
+CCTK_FILEVERSION(CactusPUGH_PUGH_Registerables_c)
+
+/********************************************************************
+ *********************  Macro Definitions  **************************
+ ********************************************************************/
+
+/********************************************************************
+ *********************  Local Data Types  ***************************
+ ********************************************************************/
+
+/********************************************************************
+ *********************  Aliased Routine Prototypes  *****************
+ ********************************************************************/
+
+/********************************************************************
+ *********************  Scheduled Routine Prototypes  ***************
+ ********************************************************************/
+
+/********************************************************************
+ *********************  Fortran Wrapper Prototypes  *****************
+ ********************************************************************/
+
+/********************************************************************
+ *********************  Local Routine Prototypes  *******************
+ ********************************************************************/
+
+static int GenerateTopologyUnavailable(int dim, 
+                                       int total_procs, 
+                                       const int *nsize, 
+                                       const int *nghostzones, 
+                                       int *nprocs);
+static void *CreateP2LTranslationUnavailable(int dim, 
+                                             int *nprocs, 
+                                             int np);
+static int DestroyP2LTranslationUnavailable(void *table);
+
+static int PhysicalToLogicalUnavailable (void *table, int proc);
+
+static int LogicalToPhysicalUnavailable(void *table, int num);
+
+/********************************************************************
+ *********************  Local Data  *********************************
+ ********************************************************************/
+
+int (*PUGHi_GenerateTopology) (int dim, 
+                               int total_procs, 
+                               const int *nsize, 
+                               const int *nghostzones, 
+                               int *nprocs) = GenerateTopologyUnavailable;
+
+void *(*PUGHi_CreateP2LTranslation) (int dim, 
+                                     int *nprocs, 
+                                     int np) = CreateP2LTranslationUnavailable;
+
+int   (*PUGHi_DestroyP2LTranslation) (void *table) = DestroyP2LTranslationUnavailable;
+int   (*PUGHi_PhysicalToLogical) (void *table, int proc) = PhysicalToLogicalUnavailable;
+int   (*PUGHi_LogicalToPhysical) (void *table, int num) = LogicalToPhysicalUnavailable;
+
+/********************************************************************
+ ************************   External Routines    ********************
+ ********************************************************************/
+
+/* Registration routines */
+
+/*@@
+   @routine    PUGH_RegisterGenerateTopology
+   @date       Sat Oct  8 19:29:25 CDT 2005
+   @author     Dylan Stark (Thomas Dramlitsch)
+   @desc 
+   Registration routine for topology generators.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+ 
+ @@*/
+int PUGH_RegisterGenerateTopology(int (*GenerateTopology) (int dim,
+                                                           int total_procs, 
+                                                           const int *nsize, 
+                                                           const int *nghostzones, 
+                                                           int *nprocs), 
+                                  const char *topologyname)
+{
+  DECLARE_CCTK_PARAMETERS;
+  
+  if (CCTK_Equals(processor_topology, topologyname))
+  {
+    if (PUGHi_GenerateTopology != GenerateTopologyUnavailable)
+    {
+      CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING,
+                 "Attempted to register topology generator %s more than once!", 
+                 topologyname);
+    }
+    else
+    {
+      PUGHi_GenerateTopology = GenerateTopology;
+      CCTK_VInfo(CCTK_THORNSTRING, "Using topology generator: %s", 
+                 topologyname);
+      
+    }
+  }
+
+  return 0;
+}
+
+/*@@
+   @routine    PUGH_RegisterP2L
+   @date       Sat Oct  8 19:29:25 CDT 2005
+   @author     Dylan Stark (Thomas Dramlitsch)
+   @desc 
+   Registration routine for physical to logical processor mappings.   
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+ 
+ @@*/
+int PUGH_RegisterP2L(void *(*CreateP2LTranslation) (int dim, 
+                                                    int *nprocs, 
+                                                    int np),
+                     int   (*DestroyP2LTranslation) (void *table),
+                     int   (*PhysicalToLogical) (void *table, 
+                                                 int proc),
+                     int   (*LogicalToPhysical) (void *table, 
+                                                 int num),
+                     const char *p2lname)
+{
+  DECLARE_CCTK_PARAMETERS;
+  
+  if (CCTK_Equals(physical2logical, p2lname))
+  {
+    if (PUGHi_CreateP2LTranslation != CreateP2LTranslationUnavailable)
+    {
+      CCTK_VWarn(0, __LINE__,__FILE__,CCTK_THORNSTRING,
+                 "Attempted to register P2L methods %s more than once!", 
+                 physical2logical);
+    }
+    else
+    {
+      PUGHi_CreateP2LTranslation  = CreateP2LTranslation;
+      PUGHi_DestroyP2LTranslation = DestroyP2LTranslation;
+      PUGHi_PhysicalToLogical     = PhysicalToLogical;
+      PUGHi_LogicalToPhysical     = LogicalToPhysical;
+
+      CCTK_VInfo(CCTK_THORNSTRING, "Using physical to logical mappings: %s",
+                 p2lname);
+    }
+  }
+ 
+  return 0;
+}
+
+/********************************************************************
+ *********************  Local Routines  *****************************
+ ********************************************************************/
+
+ /*@@
+   @routine    GenerateTopologyUnavailable
+   @date       Wed Dec 14 21:14:55 2005
+   @author     Tom Goodale
+   @desc 
+   Generates an error if the routine has not been registered.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+ 
+ @@*/
+static int GenerateTopologyUnavailable(int dim, 
+                                       int total_procs, 
+                                       const int *nsize, 
+                                       const int *nghostzones, 
+                                       int *nprocs)
+{
+  DECLARE_CCTK_PARAMETERS;
+
+  CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, 
+              "No processor topology routine registered under the name '%s'"
+              " set in the processor_topology parameter.",
+              processor_topology);
+  return 0;
+}
+
+ /*@@
+   @routine    CreateP2LTranslationUnavailable
+   @date       Wed Dec 14 21:14:55 2005
+   @author     Tom Goodale
+   @desc 
+   Generates an error if the routine has not been registered.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+ 
+ @@*/
+static void *CreateP2LTranslationUnavailable(int dim, 
+                                             int *nprocs, 
+                                             int np)
+{
+  DECLARE_CCTK_PARAMETERS;
+
+  CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, 
+              "No physical to logical mapping routine registered under the name '%s'"
+              " set in the physical2logical parameter.",
+              physical2logical);
+  return 0;
+}
+
+ /*@@
+   @routine    DestroyP2LTranslationUnavailable
+   @date       Wed Dec 14 21:14:55 2005
+   @author     Tom Goodale
+   @desc 
+   Generates an error if the routine has not been registered.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+ 
+ @@*/
+static int DestroyP2LTranslationUnavailable(void *table)
+{
+  DECLARE_CCTK_PARAMETERS;
+
+  CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, 
+              "No physical to logical mapping routine registered under the name '%s'"
+              " set in the physical2logical parameter.",
+              physical2logical);
+  return 0;
+}
+
+ /*@@
+   @routine    PhysicalToLogicalUnavailable
+   @date       Wed Dec 14 21:14:55 2005
+   @author     Tom Goodale
+   @desc 
+   Generates an error if the routine has not been registered.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+ 
+ @@*/
+static int PhysicalToLogicalUnavailable (void *table, int proc)
+{
+  DECLARE_CCTK_PARAMETERS;
+
+  CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, 
+              "No physical to logical mapping routine registered under the name '%s'"
+              " set in the physical2logical parameter.",
+              physical2logical);
+  return 0;
+}
+
+ /*@@
+   @routine    LogicalToPhysicalUnavailable
+   @date       Wed Dec 14 21:14:55 2005
+   @author     Tom Goodale
+   @desc 
+   Generates an error if the routine has not been registered.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+ 
+ @@*/
+static int LogicalToPhysicalUnavailable(void *table, int num)
+{
+  DECLARE_CCTK_PARAMETERS;
+
+  CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, 
+              "No physical to logical mapping routine registered under the name '%s'"
+              " set in the physical2logical parameter.",
+              physical2logical);
+  return 0;
+}
diff --git a/src/SetupPGV.c b/src/SetupPGV.c
index 93a2d56..149eaa0 100644
--- a/src/SetupPGV.c
+++ b/src/SetupPGV.c
@@ -63,9 +63,6 @@ static int PUGH_SetupPGExtrasSizes(int is_gf,
                                    pGExtras *this);
 static int PUGH_ComposeIJK(int dim, int *nprocs, int *pos);
 static int PUGH_DecomposeIJK(int dim, int ijk, int *nprocs, int *pos);
-static int PUGH_IntegerRoot(int number, int invpower);
-
-static int IntSort(const void *a, const void *b);
 
 /********************************************************************
  ********************* Other Routine Prototypes *********************
@@ -241,15 +238,9 @@ pConnectivity *PUGH_SetupConnectivity(int dim,
       this->neighbours[i] = this->neighbours[0]+(2*dim*i);
     }
 
-    if (CCTK_EQUALS(processor_topology, "automatic"))
-    {
-      PUGH_GenerateAlternativeTopology(dim, total_procs, nsize, nghostzones,
-                                       this->nprocs);
-    }
-    else
-    {
-      PUGH_GenerateTopology(dim, total_procs, this->nprocs);
-    }
+    PUGHi_GenerateTopology(dim, total_procs, nsize, nghostzones, 
+                             this->nprocs);
+
     PUGH_GenerateNeighbours(dim, total_procs, this->nprocs, this->neighbours,
                             this->perme);
   }
@@ -281,211 +272,6 @@ void PUGH_DestroyConnectivity(pConnectivity **conn)
 }
 
  /*@@
-   @routine    PUGH_GenerateTopology
-   @date       Fri Nov  5 11:31:21 1999
-   @author     Tom Goodale
-   @desc
-   Generate the appropriate processor topology for this processor
-   decomposition.
-   @enddesc
-   @history
-   @hdate Tue Jan 30 17:04:50 2001 @hauthor Tom Goodale
-   @hdesc Added call to integer root function and qsort to
-          avoid problems with real to integer conversions and
-          demonstrable failure of the algorithm when dealing
-          with large prime factors.
-   @endhistory
-@@*/
-int PUGH_GenerateTopology(int dim, int total_procs, int *nprocs)
-{
-  int i;
-  int used_procs;
-  int free_procs;
-  int retval;
-  int free_dims;
-
-  used_procs = 0;
-  free_procs = total_procs;
-
-  retval = 0;
-
-  free_dims = dim;
-
-  for(i=0; i < dim; i++)
-  {
-    if((nprocs[i])>0)
-    {
-      free_dims--;
-      if(used_procs)
-      {
-        used_procs *= nprocs[i];
-      }
-      else
-      {
-        used_procs = nprocs[i];
-      }
-      if (total_procs%used_procs)
-      {
-        CCTK_WARN(0, "Inconsistent PUGH topology");
-        retval = 1;
-      }
-      else
-      {
-        free_procs = total_procs/used_procs;
-      }
-    }
-  }
-
-  /* Ok calculate topology if necessary */
-  if(free_dims && ! retval)
-  {
-    /* This algorithm gives the most number of processors
-     * in the highest dimension.
-     */
-
-    int *working;
-    int root;
-    int place;
-
-    root = free_dims;
-    working = calloc(free_dims,sizeof(int));
-#ifdef DEBUG_PUGH
-    printf("Processor topology for dimension %d\n",dim);
-#endif
-
-
-    for(i = 0; i < free_dims  ; i++)
-    {
-      working[i] = PUGH_IntegerRoot(free_procs, root);
-
-      while(free_procs % working[i]) working[i]--;
-
-#ifdef DEBUG_PUGH
-      printf(" working[%d] = %d\n",i,working[i]);
-#endif
-      free_procs /= working[i];
-      root--;
-    }
-
-
-    /* The above doesn't necessarily sort them properly
-     * e.g. if one of the factors is a prime then the
-     * above will sort the 1 before the prime.
-     */
-    qsort(working,free_dims,sizeof(int),IntSort);
-
-    for(i = 0,place=0; i < dim ; i++)
-    {
-      if(nprocs[i] <= 0)
-      {
-        nprocs[i] = working[place];
-        place++;
-      }
-
-#ifdef DEBUG_PUGH
-      printf(" nprocs[%d] = %d\n",i,nprocs[i]);
-#endif
-    }
-
-    free(working);
-  }
-
-  return retval;
-}
-
-
-
- /*@@
-   @routine    PUGH_GenerateAlternativeTopology
-   @date       Thu Feb 02 17:39:21 2005
-   @author     Frank Loeffler
-   @desc
-   Generate the appropriate processor topology for this processor
-   decomposition.
-   This routine tries to decompose using the information about the actual
-   grid size. It can fail in certain situations in which it falls back
-   gracefully to the traditional decomposition giving a warning.
-   It also does not support manually set topologies and falls back in this
-   case.
-   @enddesc
-   @history
-   @endhistory
-@@*/
-int PUGH_GenerateAlternativeTopology(int dim,
-                                     int total_procs, const int *nsize,
-                                     const int *nghostzones, int *nprocs)
-{
-  int i;
-  int max_dir, max_length;
-  int free_procs = total_procs-1;
-  int used_procs = 1;
-  
-  
-  /* Nothing to decompose here */
-  if (dim == 0)
-  {
-    return PUGH_GenerateTopology(dim, total_procs, nprocs);
-  }
-  /* If there are numbers already set up or we have funny grid sizes,
-     fall back */
-  for (i = 0; i < dim; i++)
-  {
-    if (nprocs[i] || (nsize[i]<1))
-    {
-      return PUGH_GenerateTopology(dim, total_procs, nprocs);
-    }
-  }
-  /* start with a single domain */
-  for (i = 0; i < dim; i++)
-  {
-    nprocs[i] = 1;
-  }
-  /* divide as long as there are processors left */
-  while (free_procs)
-  {
-    used_procs = total_procs - free_procs;
-    /* find the direction with the longest length, which is allowed */
-    max_dir = -1;
-    max_length = 0;
-    for (i = 0; i < dim; i++)
-    {
-      /* is one part larger than the max? (then it might be a new max) */
-      if ((nsize[i] / nprocs[i] > max_length) &&
-          /* would there be at least one real point if we divide? */
-          (nsize[i] > 2*nghostzones[i]+nprocs[i]) &&
-          /* do we have enough processors left to divide in this direction? */
-          (used_procs/nprocs[i] <= free_procs))
-      {
-        max_length = nsize[i] / nprocs[i];
-        max_dir = i;
-      }
-    }
-    /* if no such direction is found: fall back giving a warning */
-    if (max_dir == -1)
-    {
-      for (i = 0; i < dim; i++)
-      {
-        nprocs[i] = 0;
-      }
-      CCTK_WARN(CCTK_WARN_COMPLAIN,
-                "Falling back to the old PUGH topology method, overwriting "
-                "eventually set manual topologies");
-      return PUGH_GenerateTopology(dim, total_procs, nprocs);
-    }
-    /* count the new direction and go on */
-    /* note: this is garanteed to decrement at least by one, since the number
-       of used processors is naturally always >= the number of processors used
-       for one dimension */
-    free_procs -= used_procs/nprocs[max_dir];
-    nprocs[max_dir]++;
-  }
-  /* success */
-  return 0;
-}
-
-
-
- /*@@
    @routine    PUGH_GenerateNeighbours
    @date       Mon Nov  8 08:15:08 1999
    @author     Tom Goodale
@@ -507,13 +293,18 @@ int PUGH_GenerateNeighbours(int dim,
   int *pos;
   int temp;
 
+  void *table;
+
+  /* Create P2L translation table. */
+  table = (void *)PUGHi_CreateP2LTranslation(dim, nprocs, total_procs);
+
   pos = malloc(dim*sizeof(int));
 
   if(dim == 0 || pos)
   {
     for(i = 0; i < total_procs; i++)
     {
-      PUGH_DecomposeIJK(dim,i,nprocs, pos);
+      PUGH_DecomposeIJK(dim, PUGHi_PhysicalToLogical(table, i), nprocs, pos);
 
       for(idim = 0; idim < dim ; idim++)
       {
@@ -522,13 +313,15 @@ int PUGH_GenerateNeighbours(int dim,
 
         if(pos[idim] > -1)
         {
-          neighbours[i][idim*2] = PUGH_ComposeIJK(dim, nprocs, pos);
+          neighbours[i][idim*2] = 
+           PUGHi_LogicalToPhysical(table, PUGH_ComposeIJK(dim, nprocs, pos));
         }
         else if(perme[idim])
         {
           temp = pos[idim];
           pos[idim] = nprocs[idim]-1;
-          neighbours[i][idim*2] = PUGH_ComposeIJK(dim, nprocs, pos);
+          neighbours[i][idim*2] = 
+            PUGHi_LogicalToPhysical(table, PUGH_ComposeIJK(dim, nprocs, pos));
           pos[idim] = temp;
         }
         else
@@ -543,13 +336,15 @@ int PUGH_GenerateNeighbours(int dim,
 
         if(pos[idim] < nprocs[idim])
         {
-          neighbours[i][idim*2+1] = PUGH_ComposeIJK(dim, nprocs, pos);
+          neighbours[i][idim*2+1] = 
+            PUGHi_LogicalToPhysical(table, PUGH_ComposeIJK(dim, nprocs, pos));
         }
         else if(perme[idim])
         {
           temp = pos[idim];
           pos[idim] = 0;
-          neighbours[i][idim*2+1] = PUGH_ComposeIJK(dim, nprocs, pos);
+          neighbours[i][idim*2+1] = 
+            PUGHi_LogicalToPhysical(table, PUGH_ComposeIJK(dim, nprocs, pos));
           pos[idim] = temp;
         }
         else
@@ -568,6 +363,9 @@ int PUGH_GenerateNeighbours(int dim,
     retval = 1;
   }
 
+  /* Destory P2L translation table. */
+  PUGHi_DestroyP2LTranslation(table);
+  
   free(pos);
 
 #ifdef DEBUG_PUGH
@@ -1593,79 +1391,3 @@ void PUGH_DestroyGArray(pGA **GA)
 /********************************************************************
  *********************     Local Routines   *************************
  ********************************************************************/
-
- /*@@
-   @routine    PUGH_IntegerRoot
-   @date       Tue Jan 30 17:06:21 2001
-   @author     Tom Goodale
-   @desc
-   Generate the highest integer below a given integer root of an integer.
-   @enddesc
-
-   @var     number
-   @vdesc   The number to take the root of
-   @vtype   int
-   @vio     in
-   @endvar
-   @var     invpower
-   @vdesc   The root to take
-   @vtype   int
-   @vio     in
-   @endvar
-
-   @returntype int
-   @returndesc
-   The highest integer below the desired root.
-   @endreturndesc
-@@*/
-static int PUGH_IntegerRoot(int number, int invpower)
-{
-  int i;
-  int tmp;
-  int root;
-
-  for(root = 1; root <= number; root++)
-  {
-    for(i=1, tmp=root; i < invpower; i++, tmp*=root);
-
-    if(tmp > number)
-    {
-      root--;
-      break;
-    }
-  }
-
-  return root;
-}
-
-
- /*@@
-   @routine    IntSort
-   @date       Tue Jan 30 17:08:47 2001
-   @author     Tom Goodale
-   @desc
-               Sorts two integers for the qsort routine.
-   @enddesc
-
-   @var        a
-   @vdesc      Pointer to first integer to compare
-   @vtype      const void *
-   @vio        in
-   @endvar
-   @var        b
-   @vdesc      Pointer to second integer to compare
-   @vtype      const void *
-   @vio        in
-   @endvar
-
-   @returntype int
-   @returndesc
-               -ve if b is greater than a.<BR>
-               +ve if a is greater than b.<BR>
-               0   if a is equal to b.
-   @endreturndesc
-@@*/
-static int IntSort (const void *a, const void *b)
-{
-  return (*(const int *) a - *(const int *) b);
-}
diff --git a/src/Topology.c b/src/Topology.c
new file mode 100644
index 0000000..457ab4e
--- /dev/null
+++ b/src/Topology.c
@@ -0,0 +1,375 @@
+ /*@@
+   @file      Topology.c
+   @date      Wed Sep 13 20:10:24 2000
+   @author    Tom Goodale
+   @desc 
+   Topology routines
+   @enddesc
+   @version  $Header$
+ @@*/
+
+#include "cctk.h"
+
+#include "pugh_Register.h"
+
+static const char *rcsid = "$Header$";
+
+CCTK_FILEVERSION(CactusPUGH_PUGH_Topology_c);
+
+/********************************************************************
+ *********************  Macro Definitions  **************************
+ ********************************************************************/
+
+/********************************************************************
+ *********************  Local Data Types  ***************************
+ ********************************************************************/
+
+/********************************************************************
+ *********************  Aliased Routine Prototypes  *****************
+ ********************************************************************/
+
+/********************************************************************
+ *********************  Scheduled Routine Prototypes  ***************
+ ********************************************************************/
+
+/********************************************************************
+ *********************  Fortran Wrapper Prototypes  *****************
+ ********************************************************************/
+
+/********************************************************************
+ *********************  Local Routine Prototypes  *******************
+ ********************************************************************/
+
+static int FranksTopology(int dim,
+                          int total_procs, 
+                          const int *nsize,
+                          const int *nghostzones, 
+                          int *nprocs);
+
+static int TraditionalTopology(int dim, 
+                               int total_procs, 
+                               const int *nsize,
+                               const int *nghostzones, 
+                               int *nprocs);
+
+static int IntegerRoot(int number, int invpower);
+
+static int IntSort (const void *a, const void *b);
+
+/********************************************************************
+ *********************  Local Data  *********************************
+ ********************************************************************/
+
+/********************************************************************
+ *********************  Aliased Routines  ***************************
+ ********************************************************************/
+
+/********************************************************************
+ *********************  Scheduled Routines  *************************
+ ********************************************************************/
+/* Fortran wrappers appear immediately after C definition of function */
+
+int PUGH_RegisterPUGHTopologyRoutines(void)
+{
+
+  PUGH_RegisterGenerateTopology(FranksTopology, 
+                                "automatic");
+  PUGH_RegisterGenerateTopology(TraditionalTopology, 
+                                "automatic_old");
+  return 0;
+}
+
+/********************************************************************
+ *********************  Other External Routines  ********************
+ ********************************************************************/
+
+/********************************************************************
+ *********************  Local Routines  *****************************
+ ********************************************************************/
+
+ /*@@
+   @routine    TraditionalTopology
+   @date       Fri Nov  5 11:31:21 1999
+   @author     Tom Goodale
+   @desc
+   Generate the appropriate processor topology for this processor
+   decomposition.
+   @enddesc
+   @history
+   @hdate Tue Jan 30 17:04:50 2001 @hauthor Tom Goodale
+   @hdesc Added call to integer root function and qsort to
+          avoid problems with real to integer conversions and
+          demonstrable failure of the algorithm when dealing
+          with large prime factors.
+   @endhistory
+@@*/
+static int TraditionalTopology(int dim, 
+                               int total_procs, 
+                               const int *nsize,
+                               const int *nghostzones, 
+                               int *nprocs)
+{
+  int i;
+  int used_procs;
+  int free_procs;
+  int retval;
+  int free_dims;
+
+  used_procs = 0;
+  free_procs = total_procs;
+
+  retval = 0;
+
+  free_dims = dim;
+
+  for(i=0; i < dim; i++)
+  {
+    if((nprocs[i])>0)
+    {
+      free_dims--;
+      if(used_procs)
+      {
+        used_procs *= nprocs[i];
+      }
+      else
+      {
+        used_procs = nprocs[i];
+      }
+      if (total_procs%used_procs)
+      {
+        CCTK_WARN(0, "Inconsistent PUGH topology");
+        retval = 1;
+      }
+      else
+      {
+        free_procs = total_procs/used_procs;
+      }
+    }
+  }
+
+  /* Ok calculate topology if necessary */
+  if(free_dims && ! retval)
+  {
+    /* This algorithm gives the most number of processors
+     * in the highest dimension.
+     */
+
+    int *working;
+    int root;
+    int place;
+
+    root = free_dims;
+    working = calloc(free_dims,sizeof(int));
+#ifdef DEBUG_PUGH
+    printf("Processor topology for dimension %d\n",dim);
+#endif
+
+
+    for(i = 0; i < free_dims  ; i++)
+    {
+      working[i] = IntegerRoot(free_procs, root);
+
+      while(free_procs % working[i]) working[i]--;
+
+#ifdef DEBUG_PUGH
+      printf(" working[%d] = %d\n",i,working[i]);
+#endif
+      free_procs /= working[i];
+      root--;
+    }
+
+
+    /* The above doesn't necessarily sort them properly
+     * e.g. if one of the factors is a prime then the
+     * above will sort the 1 before the prime.
+     */
+    qsort(working,free_dims,sizeof(int),IntSort);
+
+    for(i = 0,place=0; i < dim ; i++)
+    {
+      if(nprocs[i] <= 0)
+      {
+        nprocs[i] = working[place];
+        place++;
+      }
+
+#ifdef DEBUG_PUGH
+      printf(" nprocs[%d] = %d\n",i,nprocs[i]);
+#endif
+    }
+
+    free(working);
+  }
+
+  return retval;
+}
+
+
+
+ /*@@
+   @routine    FranksTopology
+   @date       Thu Feb 02 17:39:21 2005
+   @author     Frank Loeffler
+   @desc
+   Generate the appropriate processor topology for this processor
+   decomposition.
+   This routine tries to decompose using the information about the actual
+   grid size. It can fail in certain situations in which it falls back
+   gracefully to the traditional decomposition giving a warning.
+   It also does not support manually set topologies and falls back in this
+   case.
+   @enddesc
+   @history
+   @endhistory
+@@*/
+static int FranksTopology(int dim,
+                          int total_procs, 
+                          const int *nsize,
+                          const int *nghostzones, 
+                          int *nprocs)
+{
+  int i;
+  int max_dir, max_length;
+  int free_procs = total_procs-1;
+  int used_procs = 1;
+  
+  
+  /* Nothing to decompose here */
+  if (dim == 0)
+  {
+    return TraditionalTopology(dim, total_procs, nsize, 
+                                        nghostzones, nprocs);
+  }
+  /* If there are numbers already set up or we have funny grid sizes,
+     fall back */
+  for (i = 0; i < dim; i++)
+  {
+    if (nprocs[i] || (nsize[i]<1))
+    {
+      return TraditionalTopology(dim, total_procs, nsize, 
+                                          nghostzones, nprocs);
+    }
+  }
+  /* start with a single domain */
+  for (i = 0; i < dim; i++)
+  {
+    nprocs[i] = 1;
+  }
+  /* divide as long as there are processors left */
+  while (free_procs)
+  {
+    used_procs = total_procs - free_procs;
+    /* find the direction with the longest length, which is allowed */
+    max_dir = -1;
+    max_length = 0;
+    for (i = 0; i < dim; i++)
+    {
+      /* is one part larger than the max? (then it might be a new max) */
+      if ((nsize[i] / nprocs[i] > max_length) &&
+          /* would there be at least one real point if we divide? */
+          (nsize[i] > 2*nghostzones[i]+nprocs[i]) &&
+          /* do we have enough processors left to divide in this direction? */
+          (used_procs/nprocs[i] <= free_procs))
+      {
+        max_length = nsize[i] / nprocs[i];
+        max_dir = i;
+      }
+    }
+    /* if no such direction is found: fall back giving a warning */
+    if (max_dir == -1)
+    {
+      for (i = 0; i < dim; i++)
+      {
+        nprocs[i] = 0;
+      }
+      CCTK_WARN(CCTK_WARN_COMPLAIN,
+                "Falling back to the old PUGH topology method, overwriting "
+                "eventually set manual topologies");
+      return TraditionalTopology(dim, total_procs, nsize, 
+                                 nghostzones, nprocs);
+    }
+    /* count the new direction and go on */
+    /* note: this is garanteed to decrement at least by one, since the number
+       of used processors is naturally always >= the number of processors used
+       for one dimension */
+    free_procs -= used_procs/nprocs[max_dir];
+    nprocs[max_dir]++;
+  }
+  /* success */
+  return 0;
+}
+
+ /*@@
+   @routine    IntegerRoot
+   @date       Tue Jan 30 17:06:21 2001
+   @author     Tom Goodale
+   @desc
+   Generate the highest integer below a given integer root of an integer.
+   @enddesc
+
+   @var     number
+   @vdesc   The number to take the root of
+   @vtype   int
+   @vio     in
+   @endvar
+   @var     invpower
+   @vdesc   The root to take
+   @vtype   int
+   @vio     in
+   @endvar
+
+   @returntype int
+   @returndesc
+   The highest integer below the desired root.
+   @endreturndesc
+@@*/
+static int IntegerRoot(int number, int invpower)
+{
+  int i;
+  int tmp;
+  int root;
+
+  for(root = 1; root <= number; root++)
+  {
+    for(i=1, tmp=root; i < invpower; i++, tmp*=root);
+
+    if(tmp > number)
+    {
+      root--;
+      break;
+    }
+  }
+
+  return root;
+}
+
+ /*@@
+   @routine    IntSort
+   @date       Tue Jan 30 17:08:47 2001
+   @author     Tom Goodale
+   @desc
+               Sorts two integers for the qsort routine.
+   @enddesc
+
+   @var        a
+   @vdesc      Pointer to first integer to compare
+   @vtype      const void *
+   @vio        in
+   @endvar
+   @var        b
+   @vdesc      Pointer to second integer to compare
+   @vtype      const void *
+   @vio        in
+   @endvar
+
+   @returntype int
+   @returndesc
+               -ve if b is greater than a.<BR>
+               +ve if a is greater than b.<BR>
+               0   if a is equal to b.
+   @endreturndesc
+@@*/
+static int IntSort (const void *a, const void *b)
+{
+  return (*(const int *) a - *(const int *) b);
+}
diff --git a/src/include/pughi.h b/src/include/pughi.h
index 84cf6ba..8b4083c 100644
--- a/src/include/pughi.h
+++ b/src/include/pughi.h
@@ -16,6 +16,20 @@ extern "C"
 {
 #endif
 
+extern int (*PUGHi_DestroyP2LTranslation) (void *table);
+extern int (*PUGHi_PhysicalToLogical) (void *table, int proc);
+extern int (*PUGHi_LogicalToPhysical) (void *table, int num);
+
+extern int (*PUGHi_GenerateTopology)(int dim, 
+                                     int total_procs, 
+                                     const int *nsize,
+                                     const int *nghostzones, 
+                                     int *nprocs);
+
+extern void *(*PUGHi_CreateP2LTranslation)(int dim, 
+                                           int *nprocs, 
+                                           int np);
+
 pConnectivity *PUGH_SetupConnectivity(int dim, 
                                       int total_procs, 
                                       const int *nsize,
@@ -50,10 +64,6 @@ void PUGH_DestroyConnectivity(pConnectivity **conn);
 void PUGH_DestroyPGExtras(pGExtras **PGExtras);
 void PUGH_DestroyPGH(pGH **pughGH);
 
-int PUGH_GenerateTopology(int dim, int total_procs, int *nprocs);
-int PUGH_GenerateAlternativeTopology(int dim, int total_procs, const int *nsize,
-                                     const int *nghostzones, int *nprocs);
-
 int PUGH_GenerateNeighbours(int dim, 
                             int total_procs, 
                             int *nprocs, 
diff --git a/src/make.code.defn b/src/make.code.defn
index 3e9c6b4..4779272 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -3,4 +3,5 @@
 # Source files in this directory
 SRCS =  Startup.c GHExtension.c Overloadables.c PughUtils.c LoadAware.c       \
         Evolve.c Comm.c Storage.c SetupPGH.c SetupGroup.c SetupPGV.c          \
-        PostSendGA.c PostReceiveGA.c FinishReceiveGA.c                       
\ No newline at end of file
+        PostSendGA.c PostReceiveGA.c FinishReceiveGA.c Registration.c         \
+        P2LMappings.c Topology.c