/*
**  shm.c: 
**  a quick solution for multi-process(or) shared memory
**  programming using UNIX fork(), shmop(), semop(), and 
**  mixed Fotran and C coding. Perhaps the quickest way 
**  to convert a brute force f77 code to work in one box
**  parallelly and retain its efficiency.
**
**                                Li Ju. May 27, 1997.
*/

#include <stdio.h>
#include <errno.h>
#include <string.h>
#include <sys/ipc.h>
#include <sys/shm.h>
#include <sys/sem.h>
#include <sys/types.h>

/* define common variables */
key_t shm_KEY;   /* shared memory key */
int   shm_ID;    /* shared memory ID  */
double *shm_block;  /* shared memory logical address */
double *u,*u2,*u4,*u6,*u8,*u10,*uold,*copy;
volatile char  *schedule_flags;
/* Must be volatile because after optimization one process
   might chosen to only read schedule_flags[*My_idx-1] once into
   register, as there is no "explicit" data dependence. */
#define U_TO_U2     1
#define U2_TO_U4    2
#define U4_TO_U6    3
#define U6_TO_U8    4
#define U8_TO_U10   5
#define COLLECT_MY_DOMAIN_TO_U2  11
#define COLLECT_MY_DOMAIN_TO_U4  12
#define COLLECT_MY_DOMAIN_TO_U6  13
#define COLLECT_MY_DOMAIN_TO_U8  14
#define COLLECT_MY_DOMAIN_TO_U10 15
#define TRANSFER_FROM_V  21
#define TRANSFER_TO_V    22
#define SLAVE_DIE        99
/* addresses of passed array arguments: actually */
/* pseudonames for u's in the main program */
double **v;
double **v2;
double **v4;
double **v6;
double **v8;
double **v10;
int size_of_vector_in_bytes;
int size_of_vector_in_double;
double * address_u;

/* link with fortran common block "offset"
   to provide handles for fortran use */
extern struct
{
    int SU,SU2,SU4,SU6,SU8,SU10,SCOPY;
} offset_ ;

/* link with fortran common block "layout" */
extern struct
{
    int num_processors;
    int domain;
} layout_ ;
int *domain;

key_t sem_KEY; /* semaphore key */
int   sem_ID;  /* semaphore ID  */
volatile struct sembuf semop_add, semop_minus, semop_zero;
volatile union semun
{
    int val;
    struct semid_ds *buf;
    ushort *array;
} sem_set;

void Error(const char *fmt,...)
{
    printf("Error:%s(%s)\n",fmt,strerror(errno));
}

int generate_and_init_semaphore_ (int *num_processors)
{
    int i;
    for (sem_KEY=1;
	 ((sem_ID=semget(sem_KEY,(*num_processors),IPC_CREAT|0600))==-1)
	  &&(sem_KEY<65535);
	 sem_KEY++);
    
    if(sem_ID!=-1)
    {
	/* define semaphore operations */
	semop_add.sem_op = +1;
	semop_add.sem_flg = 0600;
	semop_minus.sem_op = -1;
	semop_minus.sem_flg = 0600;
	semop_zero.sem_op = 0;
	semop_zero.sem_flg = 0600;
	/* hypnotize the slave processes first */
	for (i=1;i<layout_.num_processors;i++)
	{
	    sem_set.val = 1;
	    semctl(sem_ID, i, SETVAL, sem_set);
	}
	return(1);
    }
    else
    {
	Error("Can not get semaphore");
	return(0);
    }
}    

void free_semaphore_ ()
{
    if(semctl(sem_ID, 0, IPC_RMID, NULL)==-1)
        printf("Unable to remove semaphore ID=%d.\n",sem_ID);
}

void slave_touch_base_ (int *My_idx)
{
    /* decreases master's semaphore */
    semop_minus.sem_num = 0;
    semop(sem_ID, &semop_minus, 1);
    /* sets his semaphore to be 1 */
    semop_add.sem_num = (ushort)(*My_idx)-1;
    semop(sem_ID, &semop_add, 1);
}

void slave_sleep_ (int *My_idx)
{
    /* wait till his semaphore is 0 */
    semop_zero.sem_num = (ushort)(*My_idx)-1;
    semop(sem_ID, &semop_zero, 1);
}

void master_set_command (char what)
{
    int i;
    /* A working session goes as:

       0) Master puts what to do in schedule_flags.

       1) Master set his semaphore to be num_processors-1,
                 set slaves' semaphores to be zero.
		 
       2) A slave wakes up when his semaphore is zero,
                 does his job defined in schedule_flags[My_idx-1],
		 decrease master's semaphore by one,
		 then sleep.

       3) Master finishes his share,
                 then sleep till his semaphore is zero. */
    
    /* put command in schedule_flags */
    memset((char *)schedule_flags,what,sizeof(char)*layout_.num_processors);
    /* set semaphore values to activate a session */
    sem_set.val = layout_.num_processors - 1;
    semctl(sem_ID, 0, SETVAL, sem_set);
    for (i=1;i<layout_.num_processors;i++)
    {
	semop_minus.sem_num = i;
	semop(sem_ID, &semop_minus, 1);
    } 
}

void master_wait ()
{
    printf ("Master wait\n");
    /* wait till his semaphore is 0 */
    semop_zero.sem_num = (ushort) 0;
    semop(sem_ID, &semop_zero, 1);
    printf ("Master finish wait\n");
}

int generate_shm_block_ (int *n3, int *num_processors)
{
    unsigned int size =
    sizeof(char)*(*num_processors)  /* memeory for schedule_flags */ 
  + sizeof(double)*2*(*n3)*(6+(*num_processors)) /* for u,..,copy */
  + sizeof(double *)*6;             /* storage for *v,*v2,*v4,*v6,*v8,*v10 */

    for (shm_KEY=1;
	 ((shm_ID=shmget(shm_KEY,size,IPC_CREAT|0600))==-1)&&(shm_KEY<65535);
	 shm_KEY++);

    if(shm_ID!=-1)
    {
	/* system will find a good logical address 
	   to attach the shared memory segment to */
	shm_block = shmat(shm_ID,(char *)0,0600);
        if(((int)shm_block)==-1)
	{
            Error("Failed to attach shared memory segment");
	    return (0);
	}
	else
	{
	    u   = (double *) shm_block;
	    u2  = u  + 2*(*n3); 
	    u4  = u2 + 2*(*n3); 
	    u6  = u4 + 2*(*n3); 
	    u8  = u6 + 2*(*n3); 
	    u10 = u8 + 2*(*n3);
	    copy = u10 + 2*(*n3);
	    v   = (double **) (copy+2*(*n3)*(*num_processors));
	    v2  = v  + 1;
	    v4  = v2 + 1;
	    v6  = v4 + 1;
	    v8  = v6 + 1;
	    v10 = v8 + 1;
	    schedule_flags = (char *)(v10+1);
	    
            /* offset list: (index in terms of double complex) */
	    offset_.SU   = 0;
	    offset_.SU2  = offset_.SU  + (*n3);
	    offset_.SU4  = offset_.SU2 + (*n3);
	    offset_.SU6  = offset_.SU4 + (*n3);
	    offset_.SU8  = offset_.SU6 + (*n3);
	    offset_.SU10 = offset_.SU8 + (*n3);
	    size_of_vector_in_double = 2*(*n3);
	    size_of_vector_in_bytes  = 2*(*n3)*sizeof(double);
	    domain = &layout_.domain;
	    /* clear the entire shared memory region */
	    memset(shm_block,0,size);
	    return (1);
	}
    }
    else
    {
	Error("Failed to get ** bytes of shared memory");
	return (0);
    }
}


/* By convention 'start' and 'finish' are
   fortran indices starting from 1, and 
   counting double precision complex numbers.
   It also corresponds to domain boundaries 
   in layout_. */

void quick_add_copy
(double *add_to_where, int start, int finish)
{
    register int i,j;
	for (j=0;j<layout_.num_processors;j++)
    for (i=(start-1)*2;i<finish*2;i++)
	    add_to_where[i] += copy[i+j*size_of_vector_in_double];
}


/* everyone does his job defined in layout_ and schedule_flags */
void do_real_job (int *My_idx)
{
    double *in, *out;
    extern void f77_sparse_multiply_ (double*,double*,int*,int*);
    extern void transfer_shmdata_from_v_and_zero_higher_order_u (int, int);
    extern void transfer_shmdata_to_v2_v10(int, int);

    /* My_idx is from 1 to num_processors */
    switch(schedule_flags[*My_idx-1])
    {  
    /* 1st phase: matrix multiplication, into my copy */
    case U_TO_U2:   in=u;  goto there;
    case U2_TO_U4:  in=u2; goto there;
    case U4_TO_U6:  in=u4; goto there;
    case U6_TO_U8:  in=u6; goto there;
    case U8_TO_U10:  in=u8; 
    there: /* my part is: */
    out = copy+(*My_idx-1)*size_of_vector_in_double;
    /* clear the work space region */
    memset (out,0,size_of_vector_in_bytes);
    /* operate in the range indicated by layout_ */
    f77_sparse_multiply_ (in,out,&domain[2*(*My_idx)-2],
			  &domain[2*(*My_idx)-1]);
    break;
    
    /* 2nd phase: after everyone is done start collecting copies */
    case COLLECT_MY_DOMAIN_TO_U2: out=u2; goto further;
    case COLLECT_MY_DOMAIN_TO_U4: out=u4; goto further;
    case COLLECT_MY_DOMAIN_TO_U6: out=u6; goto further;
    case COLLECT_MY_DOMAIN_TO_U8: out=u8; goto further;
    case COLLECT_MY_DOMAIN_TO_U10: out=u10; goto further;
    further:
    quick_add_copy (out,domain[2*(*My_idx)-2],domain[2*(*My_idx)-1]);
    break;
    
    /* At the start of master_multiply_all_: receive data from v */
    case TRANSFER_FROM_V:
	transfer_shmdata_from_v_and_zero_higher_order_u
	    (domain[2*(*My_idx)-2],domain[2*(*My_idx)-1]);
	break;
	
    /* At the end of master_multiply_all_: copy data to v2 - v10 */
    case TRANSFER_TO_V:
	transfer_shmdata_to_v2_v10
	    (domain[2*(*My_idx)-2],domain[2*(*My_idx)-1]);
	break;
    }
}


void slave_wait_for_command_ (int *My_idx)
{
    slave_sleep_ (My_idx);
    printf ("slave %d: step 1\n", *My_idx); fflush(stdout);
    while (1)
    { 
        if (schedule_flags[*My_idx-1] == SLAVE_DIE)
	{
	    slave_touch_base_(My_idx);
	    return;
	}
	else
	{
	    do_real_job (My_idx);
	    printf ("slave %d: did job %d\n", *My_idx,
		    schedule_flags[*My_idx-1]); fflush(stdout);
	    slave_touch_base_ (My_idx);
	    slave_sleep_ (My_idx);
	    printf ("slave %d: wake up with job %d\n", *My_idx,
		    schedule_flags[*My_idx-1]); fflush(stdout);
	}
    }
}


void master_multiply (char what)
{
    int i, My_idx=1;
    char command;

    /* set multiply commands for everyone */
    master_set_command (what);
    
    /* works out his own share (good master) */
    do_real_job(&My_idx);
    
    /* check how the slaves are doing */
    master_wait();

    /* if the slaves are done: collect the copies */
    switch(what)
    {
    case U_TO_U2:   command=COLLECT_MY_DOMAIN_TO_U2; goto further;
    case U2_TO_U4:  command=COLLECT_MY_DOMAIN_TO_U4; goto further;
    case U4_TO_U6:  command=COLLECT_MY_DOMAIN_TO_U6; goto further;
    case U6_TO_U8:  command=COLLECT_MY_DOMAIN_TO_U8; goto further;
    case U8_TO_U10: command=COLLECT_MY_DOMAIN_TO_U10;
    further:
    master_set_command (command);
    do_real_job(&My_idx);
    master_wait();
    }
    return;
}


/* C version of the sunroutine */
void write_shm (double *complex, int offset)
{
    /* write a complex number into shared memory */
    shm_block[(offset-1)*2] = *complex;
    shm_block[(offset-1)*2+1] = *(complex+1);
}


/* start and finish are Fortran style double complex indices */
void transfer_shmdata_from_v_and_zero_higher_order_u (int start, int finish)
{
    memcpy(u+(start-1)*2, (*v)+(start-1)*2, (finish-start+1)*sizeof(double)*2);   
    memset(u2+(start-1)*2,  0,           (finish-start+1)*sizeof(double)*2);   
    memset(u4+(start-1)*2,  0,           (finish-start+1)*sizeof(double)*2);   
    memset(u6+(start-1)*2,  0,           (finish-start+1)*sizeof(double)*2);   
    memset(u8+(start-1)*2,  0,           (finish-start+1)*sizeof(double)*2);   
    memset(u10+(start-1)*2, 0,           (finish-start+1)*sizeof(double)*2);   
}

void transfer_shmdata_to_v2_v10(int start, int finish)
{
    memcpy((*v2)+(start-1)*2,  u2+(start-1)*2,  (finish-start+1)*sizeof(double)*2);   
    memcpy((*v4)+(start-1)*2,  u4+(start-1)*2,  (finish-start+1)*sizeof(double)*2);   
    memcpy((*v6)+(start-1)*2,  u6+(start-1)*2,  (finish-start+1)*sizeof(double)*2);   
    memcpy((*v8)+(start-1)*2,  u8+(start-1)*2,  (finish-start+1)*sizeof(double)*2);   
    memcpy((*v10)+(start-1)*2, u10+(start-1)*2, (finish-start+1)*sizeof(double)*2);   
}

void master_multiply_all_
(double *vv,double *vv2,double *vv4,double *vv6,double *vv8,double *vv10,
 int *ibb,
 double *v2_ibb,double *v4_ibb,double *v6_ibb,double *v8_ibb,double *v10_ibb)
{
     int i, My_idx=1;

     /* store the Fortran array
       addresses in shared memory */
    *v   =  vv;
    *v2  =  vv2;
    *v4  =  vv4;
    *v6  =  vv6;
    *v8  =  vv8;
    *v10 =  vv10;

    /* first copy v to shared memory and zero other higher order u's */
    master_set_command(TRANSFER_FROM_V);
    printf ("master: step 1\n"); fflush(stdout);
    do_real_job(&My_idx);
    master_wait();
    printf ("master: step 2\n"); fflush(stdout);
    
    /* add in perturbation: just to indicate how to use write_shm in Fortran */
    write_shm (v2_ibb,  offset_.SU2+(*ibb));
    write_shm (v4_ibb,  offset_.SU4+(*ibb));
    write_shm (v6_ibb,  offset_.SU6+(*ibb));
    write_shm (v8_ibb,  offset_.SU8+(*ibb));
    write_shm (v10_ibb, offset_.SU10+(*ibb));
    
    /* do the multiplication on many processors */
    master_multiply (U_TO_U2);
    master_multiply (U2_TO_U4);
    master_multiply (U4_TO_U6);
    master_multiply (U6_TO_U8);
    printf ("master: step 3\n"); fflush(stdout);
    master_multiply (U8_TO_U10);
    printf ("master: step 4\n"); fflush(stdout);

    /* copy the results back to MAIN_ */
    master_set_command(TRANSFER_TO_V);
    do_real_job(&My_idx);
    master_wait();
    printf ("master: step 5\n"); fflush(stdout);
}    


void read_shm_ (double *complex, int *offset)
{ /* read in a complex number from shared memory */
    *complex = shm_block[(*offset-1)*2];
    *(complex+1) = shm_block[(*offset-1)*2+1];
}


void write_shm_ (double *complex, int *offset)
{ /* write a complex number into shared memory */
    shm_block[(*offset-1)*2] = *complex;
    shm_block[(*offset-1)*2+1] = *(complex+1);
}


void zero_shm_block_ (unsigned int *offset, int *size)
{ /* write a complex number into shared memory */
    memset(shm_block+(*offset)*2,0,(*size)*sizeof(double)*2);
}


void kill_all_slaves_ ()
{
    if (layout_.num_processors>1)
	master_set_command (SLAVE_DIE);
}


void free_shm_block_()
{
    shmctl (shm_ID,IPC_RMID,NULL);
    return;
}