#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <unistd.h>
#include <sys/types.h>
#include <sys/wait.h>
#include <sys/ipc.h>
#include <sys/shm.h>
#include <sys/sem.h>

#define N 1000
#define M N/2

typedef int (*mat_t)[N];
typedef int (*mat_blk_t)[M];

struct sembuf pop, vop;
#define P(s) semop(s,&pop,1);
#define V(s) semop(s,&vop,1);

double timespent ( struct timespec t1, struct timespec t2 )
{
   double t;

   t = 1000000000. * (double)(t2.tv_sec - t1.tv_sec);
   t += (double)(t2.tv_nsec - t1.tv_nsec);
   return t / 1000000000.;
}

mat_t createsharedmat ( int *shmid, int INIT )
{
   mat_t A;
   int i, j;

   *shmid = shmget(IPC_PRIVATE, N*N*sizeof(int), 0700|IPC_CREAT);
   A = (mat_t)shmat(*shmid,0,0);
   if (INIT) for (i=0; i<N; ++i) for (j=0; j<N; ++j) A[i][j] = rand() % 19 - 9;
   return A;
}

void matmul1 ( mat_t C, mat_t A, mat_t B )
{
   int i, j, k;

   for (i=0; i<N; ++i) {
      for (j=0; j<N; ++j) {
         C[i][j] = 0;
         for (k=0; k<N; ++k) C[i][j] += A[i][k] * B[k][j];
      }
   }
}

void blkmul2 ( mat_t C, mat_t A, mat_t B, int ctag, int shmidD, int semidD )
{
   int I, J, K, i, j, k, status, *done;
   int P[M][M];
   mat_t AA, BB, CC;
   struct timespec t1, t2;

   clock_gettime(CLOCK_REALTIME, &t1);

   I = (ctag & 4) ? 1 : 0;
   J = (ctag & 2) ? 1 : 0;
   K = (ctag & 1) ? 1 : 0;

   AA = (mat_t)(A[M*I]+M*K);
   BB = (mat_t)(B[M*K]+M*J);
   CC = (mat_t)(C[M*I]+M*J);

   for (i=0; i<M; ++i) {
      for (j=0; j<M; ++j) {
         P[i][j] = 0;
         for (k=0; k<M; ++k) P[i][j] += AA[i][k] * BB[k][j];
      }
   }

   done = (int *)shmat(shmidD, 0, 0);

   P(semidD);
   status = done[2*I+J];
   if (status == 0)
      for (i=0; i<M; ++i) for (j=0; j<M; ++j) CC[i][j] = P[i][j];
   else
      for (i=0; i<M; ++i) for (j=0; j<M; ++j) CC[i][j] += P[i][j];
   done[2*I+J] = ++status;
   V(semidD);

   shmdt(done);

   clock_gettime(CLOCK_REALTIME, &t2);
   printf("    Time taken by Process %d = %.9lf sec\n", ctag, timespent(t1,t2));

   exit(0);
}

void matmul2 ( mat_t C, mat_t A, mat_t B )
{
   int semidD[4], shmidD;
   int *done;
   int i, j, k, ctag;

   shmidD = shmget(IPC_PRIVATE, 4*sizeof(int), 0700 | IPC_CREAT);
   done = shmat(shmidD, 0, 0);
   for (i=0; i<4; ++i) done[i] = 0;
   shmdt(done);

   for (i=0; i<4; ++i) {
      semidD[i] = semget(IPC_PRIVATE, 1, 0700 | IPC_CREAT);
      semctl(semidD[i], 0, SETVAL, 1);
   }

   for (i=0; i<2; ++i) {
      for (j=0; j<2; ++j) {
         for (k=0; k<2; ++k) {
            ctag = 4*i+2*j+k;
            if (!fork())
               blkmul2(C,A,B,ctag,shmidD,semidD[2*i+j]);
         }
      }
   }

   for (i=0; i<8; ++i) wait(NULL);

   shmctl(shmidD, IPC_RMID, 0);
   for (i=0; i<4; ++i) semctl(semidD[i], 0, IPC_RMID, 0);
}

void blkmul3 ( mat_t C, mat_t A, mat_t B, int ctag, int shmidP, int shmidQ, int shmidD, int semidD )
{
   int I, J, K, i, j, k, status, *done;
   mat_t AA, BB, CC;
   mat_blk_t P, Q;
   struct timespec t1, t2;

   clock_gettime(CLOCK_REALTIME, &t1);

   I = (ctag & 4) ? 1 : 0;
   J = (ctag & 2) ? 1 : 0;
   K = (ctag & 1) ? 1 : 0;

   done = (int *)shmat(shmidD, 0, 0);
   P = (mat_blk_t)shmat(shmidP, 0, 0);
   Q = (mat_blk_t)shmat(shmidQ, 0, 0);

   AA = (mat_t)(A[M*I]+M*K);
   BB = (mat_t)(B[M*K]+M*J);
   CC = (mat_t)(C[M*I]+M*J);

   for (i=0; i<M; ++i) {
      for (j=0; j<M; ++j) {
         P[i][j] = 0;
         for (k=0; k<M; ++k)
            P[i][j] += AA[i][k] * BB[k][j];
      }
   }

   P(semidD);
   status = done[2*I+J];
   if (status == 0) done[2*I+J] = 1;
   else
      for (i=0; i<M; ++i) for (j=0; j<M; ++j)
         CC[i][j] = P[i][j] + Q[i][j];
   V(semidD);

   shmdt(done); shmdt(P); shmdt(Q);

   clock_gettime(CLOCK_REALTIME, &t2);
   printf("    Time taken by Process %d = %.9lf sec\n", ctag, timespent(t1,t2));

   exit(0);
}

void matmul3 ( mat_t C, mat_t A, mat_t B )
{
   int *done, shmidD, semidD[4];
   int shmidP[8];
   int i, j, k, ctag;

   shmidD = shmget(IPC_PRIVATE, 4*sizeof(int), 0700 | IPC_CREAT);
   done = shmat(shmidD, 0, 0);
   for (i=0; i<4; ++i) done[i] = 0;
   shmdt(done);

   for (i=0; i<8; ++i)
      shmidP[i] = shmget(IPC_PRIVATE, M*M*sizeof(int), 0700 | IPC_CREAT);

   for (i=0; i<4; ++i) {
      semidD[i] = semget(IPC_PRIVATE, 1, 0700 | IPC_CREAT);
      semctl(semidD[i], 0, SETVAL, 1);
   }

   for (i=0; i<2; ++i) {
      for (j=0; j<2; ++j) {
         for (k=0; k<2; ++k) {
            ctag = 4*i+2*j+k;
            if (!fork())
               blkmul3(C,A,B,ctag,shmidP[ctag],shmidP[ctag^1],shmidD,semidD[2*i+j]);
         }
      }
   }

   for (i=0; i<8; ++i) wait(NULL);

   shmctl(shmidD, IPC_RMID, 0);
   for (i=0; i<8; ++i) shmctl(shmidP[i], IPC_RMID, 0);
   for (i=0; i<4; ++i) semctl(semidD[i], 0, IPC_RMID, 0);
}

int mateql ( mat_t C, mat_t D )
{
   int i, j;

   for (i=0; i<N; ++i) for (j=0; j<N; ++j)
      if (C[i][j] != D[i][j]) return 0;
   return 1;
}

int main ()
{
   int shmidA, shmidB, shmidC1, shmidC2, shmidC3;
   mat_t A, B, C1, C2, C3;
   struct timespec t1, t2;

   srand((unsigned int)time(NULL));

   pop.sem_num = vop.sem_num = 0;
   pop.sem_flg = vop.sem_flg = 0;
   pop.sem_op = -1; vop.sem_op = 1;

   A = createsharedmat(&shmidA,1);
   B = createsharedmat(&shmidB,1);
   C1 = createsharedmat(&shmidC1,0);
   C2 = createsharedmat(&shmidC2,0);
   C3 = createsharedmat(&shmidC3,0);

   printf("\n+++ Sequential matrix multiplication\n");
   clock_gettime(CLOCK_REALTIME, &t1);
   matmul1(C1,A,B);
   clock_gettime(CLOCK_REALTIME, &t2);
   printf("    Time taken = %.9lf sec\n", timespent(t1,t2));

   printf("\n+++ Parallel matrix multiplication 1\n");
   clock_gettime(CLOCK_REALTIME, &t1);
   matmul2(C2,A,B);
   clock_gettime(CLOCK_REALTIME, &t2);
   printf("    Time taken = %.9lf sec\n", timespent(t1,t2));
   if (!mateql(C1,C2))
      printf("*** Parallel product does not match sequential product\n");

   printf("\n+++ Parallel matrix multiplication 2\n");
   clock_gettime(CLOCK_REALTIME, &t1);
   matmul3(C3,A,B);
   clock_gettime(CLOCK_REALTIME, &t2);
   printf("    Time taken = %.9lf sec\n", timespent(t1,t2));
   if (!mateql(C1,C3))
      printf("*** Parallel product does not match sequential product\n");

   shmdt(A); shmctl(shmidA, IPC_RMID, 0);
   shmdt(B); shmctl(shmidB, IPC_RMID, 0);
   shmdt(C1); shmctl(shmidC1, IPC_RMID, 0);
   shmdt(C2); shmctl(shmidC2, IPC_RMID, 0);
   shmdt(C3); shmctl(shmidC3, IPC_RMID, 0);

   exit(0);
}