/****   Copyright (c) 1997,1998,1999-2003  Alexandros V. Gerbessiotis
 ****
 ****   Permission to use, copy, modify, and distribute this software,
 ****   and to incorporate it, in whole or in part, into other software,
 ****   is hereby granted without fee, provided that
 ****     (1) the above copyright notice and this permission notice appear in
 ****         all copies of the source code, and the above copyright notice
 ****         appear in clearly visible form on all supporting documentation
 ****         and distribution media, and
 ****     (2) any redistribution of the software, in original
 ****         or modified form, be without fee and subject to
 ****         these same conditions.
 ****   No guarantee is offered that the code works as
 ****   advertised or is absent of any, even damaging,
 ****   side-effects; use this code on your own personal risk
 ****/
#include        "ppfai.h"
#include <stdio.h>

/***** Functions in this file
 ai2mprefix:: 2-stage parallel prefix. Operates on array inputs in general
 aitmprefix:: tree-based parallel prefix.Operates on array inputs in general
 aiprefix  :: Chooses appropriate function among ai2mprefix and aitmprefix 
              based on p,L,g. Inputs are one element arrays (variables).
 aimprefix :: Same as aiprefix but on array inputs.
 aitmmprefix:: Optimized aitmprefix but for commutative (NOT assoc)operators.
 aittmprefix::tree-based parallel prefix. A cleaneed up Ladner-Fisher variant
              probably more efficient that aitmprefix. aitmprefix is similar
              but uglier than aittmprefix.
*****/

/****  PARAMS: operator : Associative operator with prototype
                void operator (char*, char*, char*, int)
                It applies the  operator to each element
                of the second and third  argument and assigns 
                the result to the first one.
        multi    :      size of the array for parallel prefix
        from     :      source
        to       : destination (result of ppf goes here)
        nbytes   : each element of multi is of this size
****/

void ai2mprefix(void (*operator)(),int multi,char *from,char *to,int nbytes)
{
 register int   pid; 		/* processor id */
 register int   nprocs; 	/* number of processors available */
 register int   i,j,size,sze,offset; 	/* index variables */
 register int   ceil_temp,left_temp,floor_temp;
 char           *temp; 			/* buffer array */

 pid = AIPID(); nprocs = AINPROCS();

 /* split array of size multi into blocks of size either*/
 /* floor(multi/nprocs) or ceil(multi/nprocs). Assign:  */
 left_temp = LEFT(multi,nprocs);
 ceil_temp = CEIL(multi,nprocs);  /* to procs less than left_temp */
 floor_temp = FLOOR(multi,nprocs);/* to procs at least left_temp */

 /* allocate temp space and check */
 temp =  (char *) malloc(ceil_temp*nprocs*nbytes);
#ifdef SANITY_CHECK
 check_if_null((void *)temp,"ai2mprefix","temp");
#endif


 AIOREGISTER(temp,ceil_temp*nprocs*nbytes);
 AIOSYNC(temp);

 /* In first phase, each processor sends its i-th block to
  * processor i 
  */
 for (i=0;i<nprocs;i++) {
    if (i < left_temp) {
          sze= ceil_temp; /* send a ceil_temp block */
          offset=i*sze;   /* starting at... */
    } else {
          sze= floor_temp; /* send a floor_temp block */
          offset=left_temp*ceil_temp+(i-left_temp)*sze;
    }
    AIOHPPUT(i,&from[offset*nbytes],temp,pid*sze*nbytes,sze*nbytes);
 }
 AIOSYNC(temp);

 /* In second phase a local prefix operation is performed */

 /* Did pid get a ceil_temp or floor_temp block? */
 if (pid < left_temp)
       size=ceil_temp;
    else
       size=floor_temp;


 for(i=1;i<nprocs;i++) {
    operator(&temp[i*size*nbytes],&temp[(i-1)*size*nbytes],
				  &temp[i*size*nbytes],size);
 }

 /* Registration for put  operations */
 /* Register if necessary. PUSHED is defined in my_bsp.h */

 /* In third phase the computed prefices are sent back  */
 /* Where does pid send this block? */
 if (pid < left_temp) {
         sze= ceil_temp;
         offset=pid*sze;
 } else {
         sze= floor_temp;
         offset=left_temp*ceil_temp+(pid-left_temp)*sze;
 }

 for (i=0;i<nprocs;i++) {
         AIOPUT(i,&temp[i*size*nbytes],to,offset*nbytes,size*nbytes);
 }
 AIOSYNC(to);  


 /* free memory and deregister information */
 AIODEREGISTER(temp);
 /* Unregister, if necessary */
 free((void*)temp);
}


/*****  aitmprefix :: Parallel prefix using a variant
        of Ladner-Fisher's binary-tree-based algorithm on a 
        regular-tree of degree d>1.
 *****/
void aitmprefix(void (*operator)(),int multi,int degree,char *from,char *to,
                 int nbytes)
{
 register int    pid; /* processor id */
 register int    nprocs; /* number of processors */
 register int    target_pid; /* target processor in routings */
 register int    test_pid; /* indicates which processor are active */
 register int    i,j,t,t_d,shift; /* temporary variables */
 register int    size; /* size in bytes of ppf array */
 char            *temp,*tmp; /* temporary variables */
 char            *temp2;
 register int    temp_d;
 register int    base;


 pid=AIPID(); nprocs=AINPROCS();
 if ((degree< 2) || (degree>  nprocs)) degree=nprocs;
 size= nbytes * multi *sizeof(char); /* size of array to be ppf'ed */
 if (nprocs == 1 ) {
    memcpy(to,from,multi*nbytes);
    return;
 }


 /* 8 times more than needed space allocation to account for possible
    alignment problems observed in some architectures
  */
 base=size;
 temp  = (char *) malloc(degree*base);
 temp2 = (char *) malloc(degree*base);
#ifdef SANITY_CHECK
 check_if_null((void *)temp,"aitmprefix","temp");
#endif

 /* allocate storage of total size degree * size  per processor */
 /* storage is allocated this way instead of using a single     */
 /* malloc call to guarantee memory allignment                  */

 /* if source coincides with destination allocate extra memory */
 if (from == to)
   tmp = (char *) malloc(size);

 /* Check return value */
#ifdef SANITY_CHECK
 if (from == to) {
    check_if_null((void *)tmp,"aitmprefix","tmp");
 }
#endif


 /* Copy source to extra allocated space, as necessary */
 if (from == to)
    memcpy(tmp,from,size);

 /* Copy source to temp */
 memcpy(temp,from,size);

 /* Registration area */
 memcpy(temp2,from,size);
 AIOREGISTER(temp,degree*base);
 AIOSYNC(temp);


 t=degree;/* node indices differing by t/degree have same ancestor */
 t_d=1; /* t_d is always equal to t/degree */
 test_pid = pid; /* Which processors are active? */

 /* Ladner-Fisher's  upward direction stage */
 while ( t_d < nprocs) {
      target_pid = (pid/t)*t; /* Where do i route to? */
      shift= (pid-target_pid)/t_d; /* Which position in its temp */

      /* Perform first routing */
		 memcpy(temp2,temp,base*degree);
      AIOSYNC(temp);
      if ((target_pid < nprocs)  && (pid == test_pid))
        for(i=0,j=target_pid;i<degree;i++,j+=t_d)
           if (j < nprocs){
                 AIOPUT(j,temp2,temp,shift*base,size);
           }
      AIOSYNC(temp);

      /* Combine results from children, as appropriate */
      if ((target_pid < nprocs)  && (pid == test_pid)) {
        for(i=1,j=target_pid+t_d;i<degree;i++,j+=t_d)
           if (j < nprocs) {
             operator(&temp[i*base],&temp[(i-1)*base],&temp[i*base],multi);
           }
           /* If participate in next superstep, store info */
           /* into temp[0], else into to field */
           if (target_pid != pid)
             memcpy(to,&temp[(shift-1)*base],size) ;
           else
             memcpy(&temp[0],&temp[(i-1)*base],size) ;
      }

      /* Update as appropriate */
      t_d =t;
      t   *=degree;
      test_pid = target_pid;
 }


 /* Get ready for second stage */
 t   = t_d;
 t_d /= degree;
 temp_d = t_d / degree;
 target_pid = (pid/t)*t;
 shift= (pid-target_pid)/t_d;

 /* Second stage. Move downwards */

 while ( temp_d > 0 ) {
   /* Copy results to children of current node */
   if ((pid == (target_pid+shift*t_d)) && (pid != 0)) {
     for(i=1,j=pid+temp_d;i<degree;i++,j+=temp_d)
        if (j < nprocs)
           AIOPUT(j,to,temp,0,size);
   }
   AIOSYNC(temp);

   t   = t_d;
   t_d = temp_d;
   temp_d = t_d/degree;
   target_pid = (pid/t)*t;
   shift= (pid-target_pid)/t_d;

   /* Combine results */
   if ((pid == (target_pid+shift*t_d)) && (target_pid != pid) && 
       (pid >= t))
     operator(to,temp,to,multi);
 }

 /* At the leaves of the tree */
 /* act operator on them and partial results obtained thus far */
 if (pid ==0) {
   if (from == to)
     memcpy(to,tmp,size);
   else
     memcpy(to,from,size);
 } else {
   if (from == to)
     operator(to,to,tmp,multi);
    else
     operator(to,to,from,multi);
 }


 /* Unregister */
    AIODEREGISTER(temp);

 /* Deallocate memory */
 if (from == to)
   free((void*)tmp);

 free((void*)temp);
 free((void*)temp2);
}

/***** aiprefix :: Wrapper function of ai2prefix aitprefix
		   for one variable only
 *****/
void    aiprefix(void (*operator)(),char *from,char *to,int nbytes)
{
 register int    k,d;

 if (nbytes < 0) {
   AIABORT("aiprefix::size less than zero");
 } else if (nbytes ==0) {
   return;
 }

 if (__my_first_call == 0 ) {
        __my_first_call = 1;
        __my_estimate_bsp_values();
 }

 if (( __my_estimate_nprocs == 1) && (from != to)) {
   memcpy(to,from,nbytes);
   return;
 }

 k= (int) (__my_estimate_l / (__my_estimate_g));
 if ((k <= 1) || (k >= __my_estimate_nprocs)) {
   d= __my_estimate_nprocs;
  } 
  else {
    d = k;
  }
  aitmmprefix(operator,1,d,from,to,nbytes);

}

/***** aimprefix same as aiprefix but for array inputs.     
 *****/

void    aimprefix(void (*operator)(),int multi,char *from,char *to,int nbytes)
{  
 register int   k,d;

 if (nbytes < 0) {
   AIABORT("aimprefix::size less than zero");
 } else if (nbytes ==0) {
   return;
 }

 if (__my_first_call == 0 ) {
        __my_first_call = 1;
        __my_estimate_bsp_values();
 }

 if (( __my_estimate_nprocs == 1) && (from != to)) {
   memcpy(to,from,multi*nbytes);
   return;
 }

 if ((multi*nbytes) > (__my_estimate_ppf * __my_estimate_nprocs)) {
   ai2mprefix(operator,multi,from,to,nbytes);
 } 
 else {
     k= (int) (__my_estimate_l / (__my_estimate_g));
     if ((k <= 1) || (k >= __my_estimate_nprocs)) {
       d= __my_estimate_nprocs;
     } 
     else {
       d = k;
     }
   ai2mprefix(operator,multi,from,to,nbytes);
/*
     aitmmprefix(operator,multi,d,from,to,nbytes);
*/
 }
}

/* aitmmprefix. This is a butterfly based parallel prefix operator.
   Note that the operator must be commutative (not just associative)
   for the operation to work correctly.
 */

void    aitmmprefix(void (*operator)(),int multi,int degree,char *from,char *to,int nbytes)
{
 register int    pid; /* processor id */
 register int    nprocs; /* number of processors */
 register int    target_pid; /* target processor in routings */
 register int    i,j,t,t_d; /* temporary variables */
 register int    size; /* size in bytes of ppf array */
 char            *temp; /* temporary variables */
 pid=AIPID();nprocs=AINPROCS();
 if ((degree< 2) || (degree> nprocs)) degree=nprocs;

 size= nbytes * multi *sizeof(char) ; /* size of array to be ppf'ed */
 if (nprocs == 1 ) {
    memcpy(to,from,size);
    return;
 }


 /* allocate temp storage for $degree$ pointers */

 temp = (char *) malloc(degree*size);
#ifdef SANITY_CHECK
 check_if_null((void *)temp,"aitmprefix","temp");
#endif


 /* allocate storage of total size degree * size  per processor */
 /* storage is allocated this way instead of using a single     */
 /* malloc call to guarantee memory allignment                  */

 /* Copy source to temp */
 memcpy(&temp[(degree-1)*size],from,size);

 /* Registration area */
 AIOREGISTER(temp,degree*size);
 AIOSYNC(temp);

 for(i=1;i<nprocs;i*=degree) {
    AIOSYNC(temp);
    for(j=1;j<degree;j++){
       target_pid=pid+i*j;
       if (target_pid<nprocs) {
         AIOHPPUT(target_pid,&temp[(degree-1)*size],temp,
                  (degree-1-j)*size,size);
       } else {
        continue;
       }
    }
    AIOSYNC(temp);
     for(j=(degree-MIN((pid/i),(degree-1)));j<degree;j++)
        operator(&temp[j*size],&temp[(j-1)*size],&temp[j*size],multi);
  }
  memcpy(to,&temp[(degree-1)*size],size);


 /* Unregister */
    AIODEREGISTER(temp);

 /* Deallocate memory */
 free((void*)temp);
}


/****   $$ai2sprefix$$ performs parallel prefix
        in two stages. For more information read
        file brdppf.rdm

PARAMS:
        operator : Associative operator with prototype
                void operator (char*, char*, char*, int)
                It applies the  operator to each element
                of the second and third  argument and assigns 
                the result to the first one.
        multi    :      size of the array for parallel prefix
        from     :      source
        to       : destination (result of ppf goes here)
        nbytes   : each element of multi is of this size
****/

/*****EXPERIMENTAL PART NOT PRODUCTION CODE ********************/

#ifdef EXPERIMENT
void ai2sprefix(void (*operator)(),int multi,char *from,char *to,int nbytes,char *To)
{
 register int   pid; 		/* processor id */
 register int   nprocs; 	/* number of processors available */
 register int   i,j,size,sze,offset; 	/* index variables */
 register int   ceil_temp,left_temp,floor_temp;
 char           *temp; 			/* buffer array */

 pid = AIPID(); nprocs = AINPROCS();

 /* split array of size multi into blocks of size either*/
 /* floor(multi/nprocs) or ceil(multi/nprocs). Assign:  */
 left_temp = LEFT(multi,nprocs);
 ceil_temp = CEIL(multi,nprocs);  /* to procs less than left_temp */
 floor_temp = FLOOR(multi,nprocs);/* to procs at least left_temp */

 /* allocate temp space and check */
 temp =  (char *) malloc(ceil_temp*nprocs*nbytes);
#ifdef SANITY_CHECK
 check_if_null((void *)temp,"ai2mprefix","temp");
#endif


 AIOREGISTER(temp,ceil_temp*nprocs*nbytes);
 AIOSYNC(temp);

 /* In first phase, each processor sends its i-th block to
  * processor i 
  */
 for (i=0;i<nprocs;i++) {
    if (i < left_temp) {
          sze= ceil_temp; /* send a ceil_temp block */
          offset=i*sze;   /* starting at... */
    } else {
          sze= floor_temp; /* send a floor_temp block */
          offset=left_temp*ceil_temp+(i-left_temp)*sze;
    }
    AIOHPPUT(i,&from[offset*nbytes],temp,pid*sze*nbytes,sze*nbytes);
 }
 AIOSYNC(temp);

 /* In second phase a local prefix operation is performed */

 /* Did pid get a ceil_temp or floor_temp block? */
 if (pid < left_temp)
       size=ceil_temp;
    else
       size=floor_temp;


 for(i=1;i<nprocs;i++) {
    operator(&temp[i*size*nbytes],&temp[(i-1)*size*nbytes],
				  &temp[i*size*nbytes],size);
 }

 /* Registration for put  operations */
 /* Register if necessary. PUSHED is defined in my_bsp.h */
#ifdef NOTBRDPPF_PUSH
    AIOREGISTER(to,multi*nbytes);
    AIOSYNC(to);
    AIOREGISTER(To,multi*nbytes);
    AIOSYNC(To);
#endif

 /* In third phase the computed prefices are sent back  */
 /* Where does pid send this block? */
 if (pid < left_temp) {
         sze= ceil_temp;
         offset=pid*sze;
 } else {
         sze= floor_temp;
         offset=left_temp*ceil_temp+(pid-left_temp)*sze;
 }

/*
 for (i=0;i<nprocs;i++) {
         AIOPUT(i,&temp[i*size*nbytes],to,offset*nbytes,size*nbytes);
 }
 AIOSYNC(to);  
 for (i=0;i<nprocs;i++) {
         AIOPUT(i,&temp[(nprocs-1)*size*nbytes],To,offset*nbytes,size*nbytes);
 }
 AIOSYNC(To);  
*/
 for (i=0;i<nprocs;i++) {
         AIOPUT(i,&temp[i*size*nbytes],to,offset*nbytes,size*nbytes);
         AIOPUT(i,&temp[(nprocs-1)*size*nbytes],To,offset*nbytes,size*nbytes);
 }
 AIOSYNC(To);  
 AIOSYNC(to);  


 /* free memory and deregister information */
 AIODEREGISTER(temp);
 /* Unregister, if necessary */
#ifdef NOTBRDPPF_PUSH
 AIODEREGISTER(to);
#endif
 free((void*)temp);
}
#endif


/*      $aitmprefix$ performs parallel prefix on all available
        processors over the associative $operator$. Operands are
        elements of an array of dimension multi each of size $nbytes$.
        It implements a variant of Ladner-Fisher's algorithm on a
        regular-tree of degree $d>1$.

        The algorithm takes twice as much time (BSP-wise) as an one-phase
        algorithm that communicate more in each superstep (full instead
        of skewed relations.

        INTERFACE
        ---------
        As in previous function.
 */


/**** A Ladner-Fisher faithful generalization of a tree-based
      parallel prefix. It works however on ANY d-ary tree.

      Cleaned up version of aitmprefix. To keep it lean,
      it is space inefficient...
 ****/

void aittmprefix(void (*operator)(),int multi,int d,char *from,char *to,
                 int nbytes)
{
 register int    pid; /* processor id */
 register int    nprocs; /* number of processors */
 register int    i,j,k,l;
 register int    size,desp,offp,upd; /* size in bytes of ppf array */
 int             *procsize;
 char            **tmp1,**tmp2,**aux1,**aux2; /* temporary variables */


 pid=AIPID(); nprocs=AINPROCS();
 if ((d< 2) || (d>  nprocs)) d=nprocs;
 size= nbytes * multi *sizeof(char); /* size of array to be ppf'ed */
 if (nprocs == 1 ) {
    memcpy(to,from,multi*nbytes);
    return;
 }

 /*** Find log base d of nprocs. i.e find the height of the tree ***/
 upd=nprocs;l=0;
 while(upd > 1){ l++; upd=CEIL(upd,d); }
 l++; /* and the number of levels is one more than the height */
 procsize=(int *) malloc(l*sizeof(int));

 /* allocate storage of total size l*d * size  per processor */
 /* storage is allocated this way instead of using a single     */
 /* malloc call to guarantee memory allignment                  */
 tmp1 = (char **) malloc(l*sizeof(char *));
 tmp2 = (char **) malloc(l*sizeof(char *));
 aux1 = (char **) malloc(l*sizeof(char *));
 aux2 = (char **) malloc(l*sizeof(char *));
 procsize=(int *) malloc(l*sizeof(int));
 for(i=0;i<l;i++) {
   tmp1[i] = (char *) malloc(d*size);
   tmp2[i] = (char *) malloc(d*size);
   aux1[i] = (char *) malloc(size);
   aux2[i] = (char *) malloc(size);
   /* Not the best way to layout the checking of mallocs...*/
#ifdef SANITY_CHECK
   check_if_null((void *)tmp1[i],"aittmprefix","tmp1");
   check_if_null((void *)tmp2[i],"aittmprefix","tmp2");
   check_if_null((void *)aux1[i],"aittmprefix","aux1");
   check_if_null((void *)aux2[i],"aittmprefix","aux2");
#endif
 }

 memcpy(&(tmp1[0][0]),from,size);
 memcpy(&(aux1[0][0]),from,size);
 /* Ladner-Fisher's Algorithm */
 
for(j=0;j<l;j++){
   AIOREGISTER(&tmp1[j][0],d*size);
   AIOREGISTER(aux2[j],size);
}
 AIOSYNC(aux2[0]);

 procsize[0]=nprocs;
 j=1; k=nprocs; l=CEIL(k,d);

 /* First Stage : Move upwards */
 while ( k >1 ) {
   desp = pid/d;
   offp = pid%d;
/*
   AIOREGISTER(&tmp1[j][0],d*size);
*/
   AIOSYNC(&tmp1[j][0]);
   if (pid <k){
     AIOPUT(desp,aux1[j-1],&tmp1[j][0],offp*size,size);
   }
   AIOSYNC(&tmp1[j][0]);
/*
   AIODEREGISTER(&tmp1[j][0]);
*/

   if (pid < l) {
     upd=MIN(d,(k-d*pid));
     for(i=1;i<upd;i++) {
        operator(&tmp1[j][i*size],&tmp1[j][(i-1)*size],&tmp1[j][i*size],multi);
     }
     if (upd >0)
      memcpy(&aux1[j][0],&tmp1[j][(upd-1)*size],size);
     else
      memcpy(&aux1[j][0],&tmp1[j][0],size);
   }
   procsize[j]= CEIL(k,d); k=l; l=CEIL(l,d);
   j++;
 }
  /* j is the number of levels */
 j--;
  /* and now j is the height   */
 /* Nice by product:*/
 /* in pid=0, aux1[j] holds X1 o X2 o ... o Xnprocs */

 /* Second stage : Move downwards */
 while ( j > 0 ) {
   if (pid ==0)
     memcpy(aux2[j-1],tmp1[j],size);
/*
   AIOREGISTER(aux2[j-1],size);
*/
   AIOSYNC(aux2[j-1]);
   if (pid < procsize[j] ) {
     for(i=1;i<d;i++) {
       if ((pid*d+i) < procsize[j-1] ) 
         AIOPUT((pid*d+i),&tmp1[j][(i-1)*size],aux2[j-1],0,size);
      }
     if ((pid*d !=0) && (pid *d < procsize[j-1])){
         AIOPUT(pid*d,aux2[j],aux2[j-1],0,size);
     }
   }
   AIOSYNC(aux2[j-1]);
/*
   AIODEREGISTER(aux2[j-1]);
*/

   if (pid < procsize[j-1]) {
     if (j>1) upd=procsize[j-2];
       else   upd=nprocs*d;
     if (pid != 0) {
       for(i=0;i<d;i++) {
         if ( (pid*d+i) < upd )
           operator(&tmp1[j-1][i*size],aux2[j-1],&tmp1[j-1][i*size],multi);
       }
     }
     else {
      /*  !!! do nothing */
     }
   }
   
  j--; 
 }

 memcpy(to,&tmp1[0][0],size);


}
