# PD06
* Mauricio Bernuy
* Claudia Noche
* Anthony Aguilar

[ToC]

## Ejercicio 1: Ecuación de calor en una dimensión

#### En base al archivo *ec_calor_sec.cpp*

### Item a)
#### Modelo paralelizado:

<details>
<summary>Expandir Código</summary>

#include <iostream>
#include <cmath>
#include <ctime>
#include <mpi.h>

double frontera(double x, double tiempo)
{
  double limite;
  if (x < 0.5)
    limite = 100.0 + 10.0 * sin(tiempo);
  else
    limite = 75.0;
  return limite;
}

double inicial(double x, double tiempo)
{
  double limite;
  limite = 95.0;
  return limite;
}

int main(int argc, char *argv[])
{
  int my_rank, comm_sz;
  MPI_Init(NULL, NULL);
  MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);
  MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);

  int i, j, j_min = 0, j_max = 400, tag, n = 10;
  int n_com = n * comm_sz;
  double k = 0.002;
  double tiempo, dt, tmax = 10.0, tmin = 0.0, tnew;
  double u[n_com], x[n_com], dx;
  double sub_u[n], unew[n], sub_x[n];
  double x_max = 1.0, x_min = 0.0;

  dt = (tmax - tmin) / (double)(j_max - j_min);
  dx = (x_max - x_min) / (double)(n_com - 1);

  if (!my_rank)
  {
    x[0] = 0;
    for (i = 1; i < n_com; i++)
    {
      x[i] = x[i - 1] + dx;
    }
    //  Inicializacion.
    tiempo = tmin;
    for (i = 0; i < n_com; i++)
      u[i] = inicial(x[i], tiempo);
  }

  MPI_Scatter(u, n, MPI_DOUBLE, sub_u,
              n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Scatter(x, n, MPI_DOUBLE, sub_x,
              n, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  // Valores de temperatura al siguiente intervalo de tiempo
  double next_u, prev_u, next_u_recv, prev_u_recv;
  for (j = 1; j <= j_max; j++)
  {
    tnew += dt;

    if (my_rank != 0)
      MPI_Send(&sub_u[0], 1, MPI_DOUBLE, my_rank - 1, 1, MPI_COMM_WORLD);

    if (my_rank != comm_sz - 1)
      MPI_Recv(&next_u_recv, 1, MPI_DOUBLE, my_rank + 1, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE);

    int fst = 0;
    int lst = n;

    for (i = fst; i < lst; i++)
    {
      if (i < 1)
      {
        if (my_rank != 0)
        {
          MPI_Recv(&prev_u_recv, 1, MPI_DOUBLE, my_rank - 1, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
          prev_u = prev_u_recv;
        }
        else
          prev_u = 0;
      }
      else
        prev_u = sub_u[i - 1];

      if (i + 1 >= n)
      {
        if (my_rank != comm_sz - 1)
          next_u = next_u_recv;
        else
          next_u = 0;
      }
      else
        next_u = sub_u[i + 1];

      // calc unew[i]
      unew[i] = sub_u[i] + (dt * k / dx / dx) * (prev_u - 2.0 * sub_u[i] + next_u);

      // send last if needed
      if (i + 1 >= n)
      {
        if (my_rank != comm_sz - 1)

          MPI_Send(&unew[i], 1, MPI_DOUBLE, my_rank + 1, 1, MPI_COMM_WORLD);
      }
    }

    if (my_rank == 0)
      unew[0] = frontera(sub_x[0], tnew);

    if (my_rank == comm_sz - 1)
      unew[n - 1] = frontera(0, tnew);

    // Actualizar tiempo y temperatura
    tiempo = tnew;

    // update and print if last iteration
    for (i = 0; i < n; i++)
    {
      sub_u[i] = unew[i];
      if (j == j_max)
        printf("%f, %f\n", sub_x[i], sub_u[i]);
    }
  }

  MPI_Finalize();

  return 0;
}

</details>

###### > p1.cpp

#### Código en Python para generar tabla ordenada:

<details>
<summary>Expandir Código</summary>

from subprocess import Popen, PIPE
import pandas as pd
import numpy as np
import tabulate

pr = 4

print("\nUsing:", pr, "processes")
print("-------------------------------------")

df = pd.DataFrame(columns=["x","temperature"])
for pr_i in [pr]:
  program_path = f'mpirun -np {pr_i} ./a.out'
  p = Popen([program_path], shell=True, stdout=PIPE, stdin=PIPE)
  for ln in iter(p.stdout.readline,""):
    line = ln.strip().decode()
    if line == "":
      break
    x = line.split(", ")
    x = [float(i) for i in x]
    df = pd.concat([df, pd.DataFrame([x], columns=["x","temperature"])])
  

  df = df.sort_values(by=['x'])
  df.index = np.arange(1, len(df) + 1)
  print(df.to_markdown(tablefmt="grid"))

</details>

###### > tester.py

#### Tabla resultante
##### Using: 4 processes
-------------------------------------

| | x | temperature |
|:-:|:-:|:-:|
| 1 | 0 | 94.5598 |
| 2 | 0.025641 | 98.3954 |
| 3 | 0.051282 | 100.089 |
| 4 | 0.076923 | 100.376 |
| 5 | 0.102564 | 99.8944 |
| 6 | 0.128205 | 99.0992 |
| 7 | 0.153846 | 98.2676 |
| 8 | 0.179487 | 97.5369 |
| 9 | 0.205128 | 96.9532 |
| 10 | 0.230769 | 96.5122 |
| 11 | 0.25641 | 96.1879 |
| 12 | 0.282051 | 95.9413 |
| 13 | 0.307692 | 95.7555 |
| 14 | 0.333333 | 95.6115 |
| 15 | 0.358974 | 95.4973 |
| 16 | 0.384615 | 95.4061 |
| 17 | 0.410256 | 95.3344 |
| 18 | 0.435897 | 95.2811 |
| 19 | 0.461538 | 95.2456 |
| 20 | 0.487179 | 95.2281 |
| 21 | 0.512821 | 95.2287 |
| 22 | 0.538462 | 95.2458 |
| 23 | 0.564103 | 95.281 |
| 24 | 0.589744 | 95.334 |
| 25 | 0.615385 | 95.4053 |
| 26 | 0.641026 | 95.4962 |
| 27 | 0.666667 | 95.6097 |
| 28 | 0.692308 | 95.7526 |
| 29 | 0.717949 | 95.9369 |
| 30 | 0.74359 | 96.1815 |
| 31 | 0.769231 | 96.5123 |
| 32 | 0.794872 | 96.9533 |
| 33 | 0.820513 | 97.5369 |
| 34 | 0.846154 | 98.2676 |
| 35 | 0.871795 | 99.0992 |
| 36 | 0.897436 | 99.8944 |
| 37 | 0.923077 | 100.376 |
| 38 | 0.948718 | 100.089 |
| 39 | 0.974359 | 98.3954 |
| 40 | 1 | 94.5598 |

---

## Ejercicio 2:

#### En base al archivo *bubblesort.cpp*

### Item a)
#### Paralelice el código secuencial (bubblesort.cpp) con MPI u OMP (a elección) y determine Tp y speedup S

<details>
<summary>Expandir Código</summary>

#include <omp.h>

#include <iostream>
using namespace std;

#define LIMIT 1000

void imprimir_array(int A[], int N);
void Sort(int A[], int N);
void Swap(int *X, int *Y);
void leer_array(int array[], int *limite, char *filename);

int numeros[LIMIT]; // tamanho maximo del array
int size = LIMIT;

int main(int argc, char *argv[])
{
  if (argc < 2)
  { // no hay archivo => exit
    exit(1);
  }
  omp_set_num_threads(8);

  // se leen los numeros de un archivo
  leer_array(numeros, &size, argv[1]);
  // se imprime el array
  imprimir_array(numeros, size);
  // se ordena el  array
  Sort(numeros, size);
  // se imprime el array ordenado
  printf("\nArray ordenado:\n");
  imprimir_array(numeros, size);

  return 0;
}

void leer_array(int array[], int *size, char *filename)
{
  FILE *opened; // file
  char buf[30];
  int numero;      // numero de elementos leidos
  int counter = 0; // counter for elems in file

  opened = fopen(filename, "r");

  while (!feof(opened) &&
         counter < *size)
  {                            // lee hasta EOF o muchos elementos
    fgets(buf, 30, opened);    // lee una linea, max 30 caracteres
    numero = atoi(buf);        // convierte a entero
    array[counter++] = numero; // dar valores al array
    if (counter >= *size)      // muchos numeros para el array
      cout << "...muchos numeros para el array\n";
  }
  fclose(opened);

  *size = --counter; // ajusta el tamaño
  cout << *size << " numeros leidos" << endl;
}

void imprimir_array(int A[], int N)
{
  int i;
  cout << "[";
  for (i = 0; i < N; i++)
    cout << A[i] << " ";
  cout << "]\n";
}
void CompareExchange(int *X, int *Y)
{
  int t = *X ^ *Y;
  *X = t ^ max(*X, *Y); // X is minimum now
  *Y = t ^ *X;          // Y is maximum now
}

void setMaxAtEnd(int A[], int N)
{
  int m = 0;
  for (int i = 0; i < N; i++)
  {
    m = A[i] > A[m] ? i : m;
  }
  Swap(&A[m], &A[N - 1]);
}

void ParallelSort(int A[], int N)
{
  for (int i = 0; i < N; i++)
  {
#pragma omp parallel for
    for (int j = 0; j < N / 2; j++)
    {
      std::cout<< omp_get_thread_num() << "\n";

      if (i % 2)
      { // impar
        CompareExchange(&A[j * 2 + 1], &A[j * 2 + 2]);
      }
      else
      { // par
        CompareExchange(&A[j * 2], &A[j * 2 + 1]);
      }
    }
  }
}

void Sort(int A[], int N)
{
  if (N % 2)
  { // impar
    ParallelSort(A, N);
  }
  else
  {
    setMaxAtEnd(A, N);
    ParallelSort(A, N - 1);
  }
}

void Swap(int *X, int *Y)
{
  int temp = *X;
  *X = *Y;
  *Y = temp;
}

</details>

###### > p2_1.cpp

Cómo es una implementación con OMP, se considera el tiempo de comunicación como O(1), debido a la memoria compartida.

$$ T_p = O({\frac{n^2}{p}}) $$

$$ S_p = {\frac{n^2}{\frac{n^2}{p}}} = p $$

$$ E = {\frac{p}{p}} = 1 $$

### Item b)
#### Desarrolle el algoritmo de **transposición par-impar** en paralelo y determine Tp y speedup S

##### Utilizando librería: **https://github.com/swenson/sort** para la función quicksort.

<details>
<summary>Expandir Código</summary>

#include <omp.h>

#include <iostream>
#include <algorithm>

#define SORT_NAME int
#define SORT_TYPE int
#define SORT_CMP(x, y) ((x) - (y))
#include "./sort.h"

using namespace std;

#define LIMIT 1000

void imprimir_array(int A[], int N);
void ParallelSort(int A[], int N);
void Swap(int *X, int *Y);
void leer_array(int array[], int *limite, char *filename);

int numeros[LIMIT]; // tamanho maximo del array
int size = LIMIT;

int main(int argc, char *argv[])
{
  if (argc < 2)
  { // no hay archivo => exit
    exit(1);
  }
  omp_set_num_threads(8);

  // se leen los numeros de un archivo
  leer_array(numeros, &size, argv[1]);
  // se imprime el array
  imprimir_array(numeros, size);
  // se ordena el  array
  ParallelSort(numeros, size);
  // se imprime el array ordenado
  printf("\nArray ordenado:\n");
  imprimir_array(numeros, size);

  return 0;
}

void leer_array(int array[], int *size, char *filename)
{
  FILE *opened; // file
  char buf[30];
  int numero;      // numero de elementos leidos
  int counter = 0; // counter for elems in file

  opened = fopen(filename, "r");

  while (!feof(opened) &&
         counter < *size)
  {                            // lee hasta EOF o muchos elementos
    fgets(buf, 30, opened);    // lee una linea, max 30 caracteres
    numero = atoi(buf);        // convierte a entero
    array[counter++] = numero; // dar valores al array
    if (counter >= *size)      // muchos numeros para el array
      cout << "...muchos numeros para el array\n";
  }
  fclose(opened);

  *size = --counter; // ajusta el tamaño
  cout << *size << " numeros leidos" << endl;
}

void imprimir_array(int A[], int N)
{
  int i;
  cout << "[";
  for (i = 0; i < N; i++)
    cout << A[i] << " ";
  cout << "]\n";
}

void imprimir_array_range(int A[], int beg, int end)
{
  int i;
  cout << "[";
  for (i = beg; i < end; i++)
    cout << A[i] << " ";
  cout << "]\n";
}

void mergeInPlace(int A[], int start, int end, int m_start, int m_end)
{
  int L[end - start + 2];
  int R[m_end - m_start + 2];
  for (int i = start; i <= end; i++)
  {
    L[i - start] = A[i];
  }
  for (int i = m_start; i <= m_end; i++)
  {
    R[i - m_start] = A[i];
  }
  L[end - start + 1] = INT32_MAX;
  R[m_end - m_start + 1] = INT32_MAX;

  int l = 0;
  int r = 0;
  int k = start;

  while (k <= (m_end))
  {
    if (L[l] <= R[r])
    {
      A[k] = L[l];
      // cout<<L[l]<<"\n";
      l++;
    }
    else
    {
      A[k] = R[r];
      // cout<<R[r]<<"\n";
      r++;
    }
    k++;
  }
}

int getStart(int id, int p, int N)
{
  N--;
  int start = (N * id) / p;
  if (id != 0)
    start++;
  return start;
}

int getEnd(int id, int p, int N)
{
  N--;
  int end = (N * (id + 1)) / p;
  if (id == p - 1)
    end = N;
  return end;
}

void ParallelSort(int A[], int N)
{
// fase de ordenamiento individual
#pragma omp parallel
  {
    int id = omp_get_thread_num();
    int p = omp_get_num_threads();
    int start = getStart(id, p, N);
    int end = getEnd(id, p, N);

    int left_start = getStart(id + 1, p, N);
    int left_end = getEnd(id + 1, p, N);
    int right_start = getStart(id + 1, p, N);
    int right_end = getEnd(id + 1, p, N);

    int_quick_sort(&A[start], end + 1 - start);

#pragma omp barrier

    for (int i = 0; i < p; i++)
    {
      if (i % 2)
      { // impar
        if (id % 2)
        {
          if (id + 1 < p)
          {
            // printf("thread_id: %d compare left with %d\n", id, id + 1);
            mergeInPlace(A, start, end, left_start, left_end);
          }
        }
      }
      else
      { // par
        if (!(id % 2))
        {
          // printf("thread_id: %d compare right with %d\n", id, id + 1);
          mergeInPlace(A, start, end, right_start, right_end);
        }
      }
#pragma omp barrier 
   
    }
  }
 }
}

</details>

###### > p2_2.cpp

Considerando el primer término como el quicksort por proceso, y el segundo viniendo de los p merges de tamaño N/p. Al ser memoria compartida, podemos obviar el término de comunicación.

$$ T_p = O(\frac{n}{p}log(n/p)) + O(N) $$

$$ S_p = \frac{n}{1+\frac{\log \left(n/p\right)}{p}} $$

$$ E = \frac{n}{p+\log \left(n\right)-\log \left(p\right)} $$

### Item c)
#### Compare ambos resultados y analice la solución experimental con la complejidad teorica discutida en clase ¿Cuál de ellos es más escalable?

Resultado visto en clase:
$$ T_p = O(\frac{n}{p}log(n/p)) + O(N) + O(N) $$

$$ S_p = \frac{n}{2+\frac{\log \left(n/p\right)}{p}} $$

$$ E = \frac{n}{2p+\log \left(n\right)-\log \left(p\right)} $$

Como utilizamos OMP, podemos negar los tiempos de comunicación del algoritmo, lo cual nos presenta una mejora en el algoritmo de transposición par-impar (Item b).

Graficando ambas ecuaciones de Speedup con p constante (strong scaling), se puede ver que el Speedup de nuestra implementación con OMP presenta una pendiente mas favorable, por lo que presentaría mejor desempeño que el diseño original.

|  |
|:--:|
| *rojo: item b), verde: item c)*|

Por otro lado, comparando las eficiencias, a simple vista se puede notar que el denominador es mayor en la implementación que se vio en clase, por lo que se considera que la implementación presentada en esta tarea tiene mejor escalabilidad.

Cabe resaltar que los resultados del Speedup y eficiencia en el algoritmo 1 parecen ser mucho mas favorables que el diseño del transposición par-impar, pero si realizamos un ploteo de la complejidad de ambos algoritmos, podemos ver que este algoritmo depende muy fuertemente del número n de elementos, creando una parábola y elevando el tiempo de ejecución de manera desproporcionada. En cambio, la transposición par impar permite un crecimiento muy semejante al lineal, por lo cual lo seguimos consideramos como un mejor algoritmo.

|  |
|:--:|
| *verde: item a), rojo: item b)* |

---