// Davide Zanin, Alan Stocco, Alessandro Serafini e Matteo Bonotto 

#include <mpi.h>
#include <stdio.h>

float *Alloca          (int);
float *MoltMatVet      (float *, float *);
float *MoltMatPar      (float *, float *);
void   CaricaDati      (FILE  *, float *, float *, float *);
void   StampaRisultati (float *, int k); 

int n,size,rank;   // n    = classi di et… considerati
                   // size = numero di processi attivati
                   // rank = rank del processo corrente

// Funzione che alloca un array di float di dimensione dim
float *Alloca (int Dim){

 float *vet;

 vet = (float *) malloc (Dim*sizeof (float));
 
 if (vet == NULL){
  printf ("Memoria esaurita\n");
  MPI_Abort (MPI_COMM_WORLD, -1);
 }        

 return (vet);
}

// Funzione che esegue la moltiplicazione tra la marice M e l'array F
float *MoltMatVet (float *M, float *F){

 float *ris;  // Puntatore che verr… utilizzato per contenere il risultato                       
 int   i,j,r;

 // Alloca tmp
 ris = Alloca (n);

 // Esegue la moltiplicazione
 for (i = 0; i < n; i++){
  ris[i] = 0;
  for (j = 0; j < n; j++)
   ris[i] = ris[i] + F[j] * M[i*n+j];
 }
 
 
 free (F);     // F non ci serve pi— poichŠ ris Š il nuovo vettore F
 return (ris);
}

// Funzione che esegue la moltiplicazione tra due matrici in parallelo 
float  *MoltMatPar (float *M2, float *M){

 int   i,j,k,r;
 // Matrice che conterr… il risultato
 float *ris;             

 // Alloca ris 
 ris = Alloca (n*n);
 // Invia un pezzo di matrice a tutti i processi
 r = MPI_Scatter (M2, n*n/size, MPI_INT, 
                  M2, n*n/size, MPI_INT, 0, MPI_COMM_WORLD);
 // Controlla che MPI_Scatter sia stata eseguita con successo
 if (r != MPI_SUCCESS){
  printf ("Errore nell'esecuzione di MPI_Scatter\n");
  MPI_Abort (MPI_COMM_WORLD, r);  
 }    

 // Esegue la moltiplicazione in parallelo
 for (i = 0; i < n/size; i++)
  for (j = 0; j < n; j++) {
   ris[i*n+j] = 0;
   for (k = 0; k < n; k++)
    ris[i*n+j] += M2[i*n+k] * M[k*n+j];
  }

 // Ricompone la matrice
 r = MPI_Gather (ris, n*n/size, MPI_INT, 
                 ris, n*n/size, MPI_INT, 0, MPI_COMM_WORLD);
 // Controlla che MPI_Gather sia stata eseguita con successo
 if (r != MPI_SUCCESS){
  printf ("Errore nell'esecuzione di MPI_Gather\n");
  MPI_Abort (MPI_COMM_WORLD, r);  
 }

 free (M2);     // M2 non ci serve pi— poichŠ ris Š la nuova matrice M2
 return (ris);
}

// Funzione che carica i dati dal file
void CaricaDati (FILE *dati, float *F, float *M, float *M2){

 int  i;
 
 // Carica i valori di F
 for(i = 0; i < n; i++)
  fscanf (dati, " %f", &F[i]);
  
 // Carica i coefficienti di natalit…
 for(i = 0; i < n; i++)
  fscanf (dati, " %f", &M[i]);
 
 // Carica i coefficienti di mortalit…  
 for(i = 0; i < n-1; i++)
  fscanf (dati, " %f", &M[(i+1)*n+i]);           

 // Copia M in M2
 for(i = 0; i< n*n; i++)
  M2[i] = M[i];

 // Chiude il file
 fclose(dati);
}

// Funzione che stampa i risultati
void StampaRisultati (float *F, int k){
 
 int i;
 int delta;
  
 delta = 100/n;    
 printf ("  Eta'   | Femmine dopo\n         |   %d anni\n", k*delta);
 printf ("---------|-------------\n");
 for (i = 0; i < n; i++)
  printf ("%-2d - %-3d |    %-2.2f\n", delta*i, delta*(i+1), F[i]); 
 
}     
     
     
int main (int argc, char **argv){

 int   i,k,r;    // k    = numero di unit… temporali da calcolare
 float *F;       //      vettore della popolazione iniziale    
 float *M,*M2;   // M  = matrici in cui verrano sistemati i vettori p e m
                 // M2 = matrice temporanea usata per calcolare M^k
 double t0,t1;   // variabili usate per calcolare il tempo di esecuzione
 FILE   *dati;
 
 r = MPI_Init (&argc, &argv);
 if (r != MPI_SUCCESS){
  printf ("Errore di inizializzazione\n");
  MPI_Abort (MPI_COMM_WORLD, r);  
 }
 
 r = MPI_Comm_rank (MPI_COMM_WORLD, &rank);
 if (r != MPI_SUCCESS){
  printf ("Errore nell'esecuzione di MPI_Comm_rank\n");
  MPI_Abort (MPI_COMM_WORLD, r);  
 } 
 
 r = MPI_Comm_size (MPI_COMM_WORLD, &size);
 if (r != MPI_SUCCESS){
  printf ("Errore nell'esecuzione di MPI_Comm_size\n");
  MPI_Abort (MPI_COMM_WORLD, r);  
 }  
        
 t0 = MPI_Wtime (); 
 
 // Apre il file e controlla che esista
 dati = fopen ("dati","r");
 if (dati == NULL){
  if (rank == 0)        
   printf ("File 'dati' non trovato");
  MPI_Finalize ();
  return -1;
 }
 
 // Legge i valori di n e k dal file
 fscanf (dati," %d",&n);
 fscanf (dati," %d",&k);
 
 // Allocca tutti gli array e matrici necessarie
 F = Alloca (n);
 M = Alloca (n*n);
 M2 = Alloca (n*n);

 // Carica i dati dal file
 CaricaDati (dati, F, M, M2);

 // Controlla che il numero di processi attivati sia un divisore di n
 if (n%size != 0){
  if (rank == 0){
   printf ("Il numero dei processi deve essere un divisore della");
   printf (" dimensione della matrice cioŠ %d\n", n);
  } 
  MPI_Finalize ();
  return 0;
 }

 //calcolo di m^k in parallelo
 for (i = 0; i < k-1; i++)
  M2 = MoltMatPar (M2, M);

 t1 = MPI_Wtime ();

 // Calcola m^k * F e stampa i risultati
 if (rank == 0){
  F = MoltMatVet (M2, F);
  StampaRisultati (F, k);
  printf ("\nTempo di esecuzione: %f\n", t1-t0);
 }

 MPI_Finalize ();
 return 0;
}