#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <errno.h>
#include <sys/time.h>
#include <time.h>
#include <math.h>
#include "pvm3.h"
#define GRPNAME "gauss"
#define DEBUG 0

void matrix_load ( char nom[], double *tab, int me, int N, int NPROC, int tids[]) {
  FILE *f;
  int i,j,idx,idxt = 0,dest,src,msgtag = me;
	const int NBLIGNES = N/NPROC;
  double * buffer;

	if ( (buffer = malloc(N*sizeof(double))) == NULL ) {
		fprintf(stderr,"Cant malloc %d bytes\n", N*sizeof(double));
		exit (-1);
	}
  	if (me == 0){
		if ((f = fopen (nom, "r")) == NULL)
			perror ("matrix_load : fopen ");
		for (i=0; i<N; i++) {
		  for (j=0; j<N; j++) {
		    fscanf (f, "%lf", (buffer+j) );
		  } 
			idx = i % NPROC;
			if (idx == 0) {
				memcpy(tab+(idxt++)*N, buffer, N*sizeof(double));
			} else {
				pvm_initsend( PvmDataDefault );
				pvm_pkdouble( buffer, N, 1);
				pvm_send( tids[idx], i );
			}
		}
		fclose (f);
	} else {
		for (i=0; i<NBLIGNES; i++, msgtag += NPROC) {
			pvm_recv( tids[0], msgtag );
      			pvm_upkdouble( buffer, N, 1 );
			idx = i % NPROC;
			memcpy(tab+(idxt++)*N, buffer, N*sizeof(double));
		}
	}
}

void matrix_save ( char nom[], double *tab, int me, int N, int NPROC, int tids[] ) {
  FILE *f;
  int i, j, idx, idxt = 0, src, dest, msgtag = me;
	const int NBLIGNES = N/NPROC;
	double *buffer;
	if ( (buffer = malloc(N*sizeof(double))) == NULL ) {
		fprintf(stderr,"Cant malloc %d bytes\n", N*sizeof(double));
		exit (-1);
	}

	if (me == 0) {
		if ((f = fopen (nom, "w")) == NULL)
			perror ("matrix_save : fopen ");
		for (i=0; i<N; i++) {
			idx = i%NPROC;
			if (idx == 0) {
				memcpy(buffer,tab+(idxt++)*N,N*sizeof(double));
			} else {
				pvm_recv( tids[idx], i );
				pvm_upkdouble( buffer, N, 1 );
			}
		for (j=0; j<N; j++) {
	    		fprintf (f, "%8.2f\t", buffer[j] );
	  	}
  		fprintf (f, "\n");
  	}
  	fclose (f);
	} else {
		for (i=0; i<NBLIGNES; i++, msgtag += NPROC) {
			pvm_initsend( PvmDataDefault );
			pvm_pkdouble( tab+i*N, N, 1);
			pvm_send( tids[0], msgtag );
		}
	}
}

void matrix_display ( double *tab, int  N, int NPROC) {
	int i,j;
	const int NBLIGNES = N/NPROC;
	for (i=0; i<NBLIGNES*N; i++) {
		printf ("%8.2f ", tab[i] );
	printf ("\n");
	}
}

void gauss (int me, double * tab, int N , int NPROC, int tids[]) {
	int i,j,k,l;
	double pivot;
	double buffer[N];
	int decalage = 0;

	for(k=0;k<(N-1);k++) {
		if(k%NPROC==me) {
			memcpy(&buffer[k], tab + (k/NPROC)*N + k,(N-k)*sizeof(double));
			pvm_initsend(PvmDataDefault);
			pvm_pkdouble(&buffer[k],N-k,1);
			for(l=0; l<NPROC;l++){
				if(l != me){
					pvm_send(tids[l], k);
				}
			}
		} else {
			pvm_recv(tids[k%NPROC], k);
			pvm_upkdouble(&buffer[k],N-k,1);
		}
		decalage = (k/NPROC) + ((me>(k%NPROC))?0:1);
		for(i=decalage;i<(N/NPROC);i++) {
			if ( fabs(buffer[k]) <= 1.0e-11 ) {
      				printf ("ATTENTION: pivot %d presque nul: %g\n", k, buffer[k] );
      				exit (-1);
			}
			pivot = (*(tab + i * N + k)/buffer[k]);
			for(j=k;j<N;j++){
				*(tab + i * N +j) = *(tab + i * N + j) - pivot*buffer[j];
			}
		}
	}
}

main(int argc, char ** argv) {
		if(DEBUG) printf("JE COMMENCE !\n");
    int NPROC = 8;              /* default nb of proc */
    int mytid;                  /* my task id */
    int *tids;                  /* array of task id */
    int me;                     /* my process number */
    int i;
    /* enroll in pvm */
    mytid = pvm_mytid();

    /* determine the size of my sibling list */
    NPROC = pvm_siblings(&tids); 
		if (DEBUG) printf("NPROC: %d\n", NPROC);
    /* WARNING: tids are in order of spawning, which is different from
       the task index JOINING the group */

    me = pvm_joingroup( GRPNAME ); /* me: task index in the group */
	if (DEBUG) fprintf(stderr,"me: %d\n", me);
    pvm_barrier( GRPNAME, NPROC );
    pvm_freezegroup ( GRPNAME, NPROC );
    for ( i = 0; i < NPROC; i++) tids[i] = pvm_gettid ( GRPNAME, i); 

/*--------------------------------------------------------------------------*/
/*           all the tasks are equivalent at that point                     */

     dowork(argc, argv, me, tids, NPROC );

     pvm_lvgroup( GRPNAME );
     pvm_exit();
}

/* Simple example passes a token around a ring */

dowork(int argc, char ** argv, int me, int tids[], int NPROC ) {
     
  int N, i, j, k;
  double *tab, pivot;
  char nom[127], nomf[127];
  FILE *f;
  struct timeval tv1, tv2;      /* for timing */
  int taille_block;
  double duree;
  
  
  if (argc != 3){
    printf ("Usage: %s <matrix size> <matrix name>\n", argv[0]);
    exit (-1);
	}

  N = atoi ( argv[1] );
  strcpy(nom, argv[2]);
  if ( (tab=malloc(N*N/NPROC*sizeof(double))) == NULL ) {
    printf ("Cant malloc %d bytes\n", N*N*sizeof(double));
    exit (-1);
  }

  if (me == 0) {
    strcpy(nomf, "pvm3/result");
    sprintf ( nomf+strlen(nomf), "%d", N );
    sprintf ( nomf+strlen(nomf), ".txt" );
    if((f = fopen (nomf, "a")) == NULL) perror ("matrix_save : fopen ");
  }
  
  gettimeofday( &tv1, (struct timezone*)0 );
  if(DEBUG) printf("Before load: me %d\n",me);
  matrix_load ( nom, tab, me, N, NPROC, tids);
  if(DEBUG) printf("After load: me %d\n",me);
  gettimeofday( &tv2, (struct timezone*)0 );
  duree = (tv2.tv_sec - tv1.tv_sec) * 1000000.0 + tv2.tv_usec - tv1.tv_usec;
  //if(me == 0) fprintf (f, "nproc %d\tloading time: %10.8f msec.\n", NPROC, duree/1000.0 );

  gettimeofday( &tv1, (struct timezone*)0 );
  gauss ( me, tab, N, NPROC, tids);
  gettimeofday( &tv2, (struct timezone*)0 );
  duree = (tv2.tv_sec - tv1.tv_sec) * 1000000 + tv2.tv_usec - tv1.tv_usec;
  if (me == 0) {
	fprintf (f, "nproc %d\tcommuting time: %10.8f msec.\n", NPROC, duree/1000.0 );
    fclose(f);
  }
	
  sprintf ( nom+strlen(nom), ".result" );
  gettimeofday( &tv1, (struct timezone*)0 );
  if(DEBUG) printf("Before save: me %d\n",me);
  matrix_save ( nom, tab, me, N, NPROC, tids);
  if(DEBUG) printf("After save: me %d\n",me);
  gettimeofday( &tv2, (struct timezone*)0 );
  
}