%%bash
lscpu

%%bash
sudo lshw -C memory

%%bash
uname -ar #r for kernel, a for all

"""
Explanation: (PERFBLAS)=
5.2 Herramientas de lenguajes de programación y del sistema operativo para perfilamiento e implementaciones de BLAS
```{admonition} Notas para contenedor de docker:
Comando de docker para ejecución de la nota de forma local:
nota: cambiar &lt;ruta a mi directorio&gt; por la ruta de directorio que se desea mapear a /datos dentro del contenedor de docker y &lt;versión imagen de docker&gt; por la versión más actualizada que se presenta en la documentación.
docker run --rm -v &lt;ruta a mi directorio&gt;:/datos --name jupyterlab_optimizacion_2 -p 8888:8888 -d palmoreck/jupyterlab_optimizacion_2:&lt;versión imagen de docker&gt;
password para jupyterlab: qwerty
Detener el contenedor de docker:
docker stop jupyterlab_optimizacion_2
Documentación de la imagen de docker palmoreck/jupyterlab_optimizacion_2:&lt;versión imagen de docker&gt; en liga.
```

```{admonition} Al final de esta nota la comunidad lectora:
:class: tip


Aprenderá diferentes herramientas para perfilamiento de código y dará interpretación a sus resultados.


Conocerá algunas operaciones de BLAS en los tres niveles. En específico se utilizará OpenBLAS para NumPy y R.


```
En esta nota revisamos algunas herramientas de los lenguajes de programación y de los sistemas GNU/Linux para perfilamiento de código. También se revisa la implementación OpenBLAS de la especificación de Basic Linear Algebra Subprograms (BLAS) en sus diferentes niveles con ejemplos de paquetes en los lenguajes.
```{margin}
Utilizamos la máquina m4.16xlarge pues en AWS no todas las instancias se pueden obtener algunas métricas con perf. Ver event-cache-misses-not-supported-by-perf-in-aws, ubuntu-12-10-perf-stat-not-supported-cycles, the-pmcs-of-ec2
```
Se presentan códigos y sus ejecuciones en una máquina m4.16xlarge con una AMI ubuntu 20.04 - ami-042e8287309f5df03 de la nube de AWS. Se utilizó en la sección de User data el script_profiling_and_BLAS.sh
La máquina m4.16xlarge tiene las siguientes características:
End of explanation
"""


import math
import time

import numpy as np
from pytest import approx
from scipy.integrate import quad

def Rcf(f,a,b,n):
    """
    Compute numerical approximation using rectangle or mid-point
    method in an interval.
    Nodes are generated via formula: x_i = a+(i+1/2)h_hat for
    i=0,1,...,n-1 and h_hat=(b-a)/n
    Args:
    
        f (float): function expression of integrand.
        
        a (float): left point of interval.
        
        b (float): right point of interval.
        
        n (int): number of subintervals.
        
    Returns:
    
        sum_res (float): numerical approximation to integral
            of f in the interval a,b
    """
    h_hat = (b-a)/n
    nodes = [a+(i+1/2)*h_hat for i in range(n)]
    sum_res = 0
    for node in nodes:
        sum_res = sum_res+f(node)
    return h_hat*sum_res


"""
Explanation: ```{admonition} Observación
:class: tip
En la celda anterior se utilizó el comando de magic %%bash. Algunos comandos de magic los podemos utilizar también con import. Ver ipython-magics
```
```{admonition} Comentario
Lo más recomendable para realizar perfilamiento de códigos es realizarlos no en jupyter notebooks. A continuación se realizarán algunas mediciones en este notebook para propósitos didácticos pero se sugiere perfilar en el intérprete de Python o ejecutando un script.py.
```
Perfilamiento: medición de tiempo en Python y IPython
Lo más natural que podemos pensar en medir es el tiempo de ejecución de nuestros códigos. Python y IPython tienen herramientas para este propósito.
Módulo: time
Imports y funciones de apoyo:
End of explanation
"""


f = lambda x: math.exp(-x**2) #using math library
a = 0
b = 1
obj, err = quad(f, a, b)
print(obj)

"""
Explanation: Objetivo:
End of explanation
"""


n = 10**6
start_time = time.time()
res = Rcf(f,a,b,n)
end_time = time.time()
secs = end_time-start_time
print("Rcf tomó",secs,"segundos" )

"""
Explanation: Medición de tiempo:
End of explanation
"""


print(res == approx(obj))

"""
Explanation: Prueba que se resuelve correctamente el problema:
End of explanation
"""


%time Rcf(f,a, b,n)

"""
Explanation: Comando de magic: %time
Nos regresa las mediciones siguientes:


CPU times que contiene:


user: mide la cantidad de tiempo de los statements que la CPU gastó para funciones que no están relacionadas con el kernel del sistema, por ejemplo ejecución de instrucciones que involucran aritmética.


sys: mide la cantidad de tiempo de los statements que la CPU gastó en funciones a nivel de kernel del sistema, por ejemplo alojamiento, lectura y escritura de variables en memoria, las relacionadas con el I/O de disco o network.


total: suma entre el user y sys para todos los cores. 




Wall time: mide el wall clock o elapsed time que se refiere al tiempo desde que inicia la ejecución de los statements hasta su finalización.


Out: resultado.
End of explanation
"""


%timeit -n 5 -r 10 Rcf(f,a, b,n)

"""
Explanation: ```{admonition} Comentarios


Para mediciones de tiempos que involucran cómputo en paralelo es posible que total exceda a Wall time por una diferencia notable pues es la suma de tiempos para todos los cores.


La diferencia entre total y Wall time da una idea de la cantidad de tiempo que se ocupó el sistema en la ejecución de statements que no involucran user ni a sys (por ejemplo un statement del tipo sleep) o en tareas no relacionadas con los statements medidos (por ejemplo si en el momento de la medición se estaban corriendo otros procesos).


```
Al medir tiempos de ejecución, siempre hay variación en la medición por lo que se recomienda realizar las mediciones un número repetido de veces. Para este caso Timeit nos ayuda. 
Timeit
Se ejecuta desde la línea de comandos, con el comando de magic %timeit o realizando import.
A continuación se mide el tiempo de ejecución para la función Rcf. Se promedian los tiempos de $n=5$ ejecuciones y se calcula su desviación estándar. Se repite lo anterior $r=10$ veces y se reporta el mejor resultado.
End of explanation
"""


import cProfile

n=10**6
cprof = cProfile.Profile()
cprof.enable()
res = Rcf(f,a,b,n)
cprof.disable()
cprof.print_stats(sort='cumtime')

"""
Explanation: ```{admonition} Observación
:class: tip
$ms$ es milisecond, $\mu s$ es microsecond y $ns$ es nanosecond.
```
```{admonition} Comentarios


%timeit se recomienda usar para secciones de código pequeñas. Si se observa una variación en las mediciones de tiempo entre distintas repeticiones entonces hay que realizar más repeticiones hasta tener un resultado estable.


%timeit desabilita temporalmente el garbage collector de Python (esto es, no habrá desalojamiento en memoria de objetos de Python que no se utilicen). Si el garbage collector es invocado en tus statements, esto puede ser una razón de posibles diferencias que se obtengan en las mediciones de tiempo.


Para su uso con un import ver timeit y how-to-use-timeit-module.


```
```{admonition} Ejercicio
:class: tip
Implementar la regla de Simpson compuesta para aproximar la misma integral y medir los tiempos de ejecución con un mismo número de subintervalos vía Timeit utilizando import. Revisar que se está aproximando correctamente el valor numérico con quad.
```
cProfile
cProfile está en la standard-library de Python como built-in. Se utiliza con la implementación CPython de Python para medir el tiempo de ejecución de cada función en el programa. Se ejecuta desde la línea de comandos, con un comando de magic o realizando import.
End of explanation
"""


cprof.dump_stats("Rcf_stats")

import pstats

p_rcf_stats = pstats.Stats("Rcf_stats")
print(p_rcf_stats.sort_stats("cumulative").print_stats(10))

"""
Explanation: El output de cProfile muestra:


El tiempo total de ejecución, el cual incluye el tiempo del bloque de código que estamos midiendo y el overhead al usar cProfile. Por esta razón se tiene un mayor tiempo de ejecución que con las mediciones de tiempo anteriores.


La columna ncalls que como el nombre indica, muestra el número de veces que se llamó a cada función. En este caso las funciones lambda y math.exp son las que se llaman un mayor número de veces: $n=10^6$ veces. 


La columna tottime muestra el tiempo que tardaron estas funciones en ejecutarse (sin llamar a otras funciones).


La columna percall es el cociente entre tottime y ncalls.


La columna cumtime contiene el tiempo gastado en la función y en las demás que llama. Por ejemplo la función Rcf llama a listcomp por lo que es natural que Rcf esté más arriba en el output ordenado de cProfile. Esto también ocurre con lambda y math.exp pues la primera llama a la segunda.


La columna de percall es un cociente entre la columna cumtime y el conteo del número de veces que se llamaron a funciones primitivas o también nombradas built in functions.


La última columna indica información de la función y la línea en la que se encuentra dentro del código. Por ejemplo para el caso de listcomp sale la línea 23 y si se va a la celda en donde se definió Rcf corresponde tal línea al llamado nodes=[a+(i+1/2)*h_hat for i in range(n)]. 


````{admonition} Comentarios


El comando de magic para ejecutar cProfile en una celda del jupyter notebook es: %prun -s cumulative Rcf(f,a, b,n)


cProfile puede ejecutarse desde la línea de comando con: python3 -m cProfile -s cumulative Rcf.py creando Rcf.py por ejemplo con:


```python
%%file Rcf.py
import math
from pytest import approx
from scipy.integrate import quad
def Rcf(f,a,b,n):
    """
    Compute numerical approximation using rectangle or mid-point
    method in an interval.
    Nodes are generated via formula: x_i = a+(i+1/2)h_hat for
    i=0,1,...,n-1 and h_hat=(b-a)/n
    Args:
    f (float): function expression of integrand.

    a (float): left point of interval.

    b (float): right point of interval.

    n (int): number of subintervals.

Returns:

    sum_res (float): numerical approximation to integral
        of f in the interval a,b
"""
h_hat = (b-a)/n
nodes = [a+(i+1/2)*h_hat for i in range(n)]
sum_res = 0
for node in nodes:
    sum_res = sum_res+f(node)
return h_hat*sum_res

if name == "main":
    a = 0
    b = 1
    f = lambda x: math.exp(-x2)
    obj, err = quad(f, a, b)
    print("objetivo: {:0.6e}".format(obj))
    n = 106
    res = Rcf(f,a,b,n)
    print("aproximación: {:0.6e}".format(res))
    print(res == approx(obj))
```
La flag -s con la ejecución vía línea de comando, indica que se ordene el resultado por el tiempo acumulado dentro de cada función.
````
El output de cProfile no es muy flexible por lo que tenemos el módulo pstats con la clase Stats para manipular los resultados dentro de Python. Por ejemplo para ordenar por cumulative time e imprimir sólo las 10 líneas más significativas se puede realizar:
```{margin}
Creamos el archivo Rcf_stats para ser usado por pstats.
```
End of explanation
"""


print(p_rcf_stats.sort_stats("cumulative").print_stats("lambda|listcomp|math"))

"""
Explanation: Si queremos filtrar los resultados para los que hagan match con las palabras lambda o listcomp o math tomadas como expresiones regulares en la columna filename:lineno(function):
End of explanation
"""


print(p_rcf_stats.strip_dirs().sort_stats("cumulative").print_stats("lambda|listcomp|math"))

"""
Explanation: Y podemos usar strip_dirs para eliminar el nombre de la ruta para mayor legibilidad:
End of explanation
"""


print(p_rcf_stats.prim_calls)

"""
Explanation: El número de veces que se llamaron a funciones primitivas o built in es:
End of explanation
"""


p_rcf_stats.strip_dirs().sort_stats("cumulative").print_callers()

p_rcf_stats.strip_dirs().sort_stats("cumulative").print_callees()

"""
Explanation: También pstats nos permite obtener qué funciones llaman a otras o son llamadas:
End of explanation
"""


p_rcf_stats.strip_dirs().sort_stats("cumulative").print_callers(10)

p_rcf_stats.strip_dirs().sort_stats("cumulative").print_callees("Rcf|lambda")

"""
Explanation: Análogamente podemos establecer los filtros usados en print_stats:
End of explanation
"""


def myfunc():
    A = np.array([[1, 2, 0],
                  [3, 4, 5]], dtype=float)
    a_1 = A[:, 0]
    return a_1

%prun -s cumulative myfunc()

"""
Explanation: ```{admonition} Comentario
cProfile sólo mide tiempos para llamadas a funciones de forma explícita no métodos que son llamados de forma interna como en el siguiente ejemplo en el que se hace una extracción de la primer columna de un NumPy array.
```
End of explanation
"""


import line_profiler

line_prof = line_profiler.LineProfiler()
n = 10**6
print(line_prof(Rcf)(f,a,b,n))

print(line_prof.print_stats())

"""
Explanation: ```{admonition} Comentario
cProfile puede usarse en compañía del paquete SnakeViz para visualizar los resultados.
```
```{admonition} Ejercicio
:class: tip
Perfilar la regla de Simpson compuesta que se implementó en el primer ejercicio con cProfile utilizando su import y  el paquete pstats en una máquina de AWS con las mismas características que la que se presenta en esta nota.
```
line_profiler
Nos ayuda a perfilar el código línea por línea. 
Una buena idea es perfilar el código primero con cProfile para identificar aquellas funciones que gastan un mayor tiempo de ejecución y posteriormente perfilarlas con line_profiler.
Se ejecuta desde la línea de comandos, con un comando de magic o realizando import.
End of explanation
"""


def Rcf_2(f,a,b,n):
    """
    Compute numerical approximation using rectangle or mid-point
    method in an interval.
    Nodes are generated via formula: x_i = a+(i+1/2)h_hat for
    i=0,1,...,n-1 and h_hat=(b-a)/n
    Args:
    
        f (float): function expression of integrand.
        
        a (float): left point of interval.
        
        b (float): right point of interval.
        
        n (int): number of subintervals.
        
    Returns:
    
        sum_res (float): numerical approximation to integral
            of f in the interval a,b
    """
    h_hat = (b-a)/n
    nodes = (a+(i+1/2)*h_hat for i in range(n))
    sum_res = 0
    for node in nodes:
        sum_res = sum_res+f(node)
    return h_hat*sum_res


line_prof.add_function(Rcf_2)
print(line_prof(Rcf_2)(f, a, b, n))

"""
Explanation: El output de line_profiler muestra:


Line #: el número de la línea.


Hits: el número de veces que la línea fue ejecutada.


Time el tiempo que tomó la ejecución de la línea. Las unidades en el que se mide el tiempo está en la primer línea: Timer unit.


Per Hit: el tiempo promedio que toma la ejecución de la línea una vez en timer's units.


\% Time: el porcentaje de tiempo que toma la ejecución de la línea correspondiente al tiempo total de ejecución de la función.


Line contents: el código a perfilar.


```{admonition} Comentarios


Para cambiar Timer unit puede usarse: line_prof.print_stats(output_unit=1e-03).


Con add_function de line_profiler es posible perfilar varias funciones.


```
```{margin}
Ver glossary: generator para la definición de un generator.
```
Utilizando generators:
End of explanation
"""


f_np = lambda x: np.exp(-x**2)

def Rcf_3(f,a,b,n):
    """
    Compute numerical approximation using rectangle or mid-point
    method in an interval.
    Nodes are generated via formula: x_i = a+(i+1/2)h_hat for
    i=0,1,...,n-1 and h_hat=(b-a)/n
    Args:
    
        f (float): function expression of integrand.
        
        a (float): left point of interval.
        
        b (float): right point of interval.
        
        n (int): number of subintervals.
        
    Returns:
    
        sum_res (float): numerical approximation to integral
            of f in the interval a,b
    """
    h_hat = (b-a)/n
    aux_vec = np.linspace(a, b, n+1)
    nodes = (aux_vec[:-1]+aux_vec[1:])/2
    return h_hat*np.sum(f(nodes))


line_prof.add_function(Rcf_3)
print(line_prof(Rcf_3)(f_np, a, b, n))

print(line_prof.print_stats())

"""
Explanation: Utilizando funcionalidad de NumPy:
End of explanation
"""


%%file Rcf_line_profiler.py
import math

from pytest import approx
from scipy.integrate import quad

@profile #decorate the functions you want to profile with kernprof 
def Rcf(f,a,b,n):
    """
    Compute numerical approximation using rectangle or mid-point
    method in an interval.
    Nodes are generated via formula: x_i = a+(i+1/2)h_hat for
    i=0,1,...,n-1 and h_hat=(b-a)/n
    Args:
    
        f (float): function expression of integrand.
        
        a (float): left point of interval.
        
        b (float): right point of interval.
        
        n (int): number of subintervals.
        
    Returns:
    
        sum_res (float): numerical approximation to integral
            of f in the interval a,b
    """
    h_hat = (b-a)/n
    nodes = [a+(i+1/2)*h_hat for i in range(n)]
    sum_res = 0
    for node in nodes:
        sum_res = sum_res+f(node)
    return h_hat*sum_res

if __name__ == "__main__":
    a = 0
    b = 1
    f = lambda x: math.exp(-x**2)
    obj, err = quad(f, a, b)
    print("objetivo: {:0.6e}".format(obj))
    n = 10**6
    res=Rcf(f,a,b,n)
    print("aproximación: {:0.6e}".format(res))
    print(res == approx(obj))

%%bash
$HOME/.local/bin/kernprof -l -v Rcf_line_profiler.py

"""
Explanation: ```{admonition} Comentarios


Para su uso con el comando de magic hay que cargarlo como una extensión de jupyter: %load_ext line_profiler y posteriormente usar: %lprun -f Rcf Rcf(f, a, b, n).


Para su uso vía la línea de comandos se debe ejecutar: kernprof -v -l &lt;script.py&gt; o bien python3 &lt;script.py&gt; como se muestra a continuación. La opción -v es view the results, -l indica que se haga el perfilamiento línea a línea. Por default se escribe un archivo con extensión .lprof el cual puede revisarse con python3 -m line_profiler &lt;script.py.lprof&gt;


```
End of explanation
"""


%%file Rcf_line_profiler_2.py
import math

import line_profiler
from pytest import approx
from scipy.integrate import quad

def Rcf(f,a,b,n):
    """
    Compute numerical approximation using rectangle or mid-point
    method in an interval.
    Nodes are generated via formula: x_i = a+(i+1/2)h_hat for
    i=0,1,...,n-1 and h_hat=(b-a)/n
    Args:
    
        f (float): function expression of integrand.
        
        a (float): left point of interval.
        
        b (float): right point of interval.
        
        n (int): number of subintervals.
        
    Returns:
    
        sum_res (float): numerical approximation to integral
            of f in the interval a,b
    """
    h_hat = (b-a)/n
    nodes = [a+(i+1/2)*h_hat for i in range(n)]
    sum_res = 0
    for node in nodes:
        sum_res = sum_res+f(node)
    return h_hat*sum_res

if __name__ == "__main__":
    a = 0
    b = 1
    f = lambda x: math.exp(-x**2)
    obj, err = quad(f, a, b)
    print("objetivo: {:0.6e}".format(obj))
    n = 10**6
    res = Rcf(f,a,b,n)
    print("aproximación: {:0.6e}".format(res))
    print(res == approx(obj))
    line_prof = line_profiler.LineProfiler()
    print(line_prof(Rcf)(f,a,b,n))
    line_prof.print_stats()    

%%bash
python3 Rcf_line_profiler_2.py

"""
Explanation: Utilizando python3 &lt;script.py&gt;:
End of explanation
"""


from memory_profiler import memory_usage

t = (Rcf, (f, a, b, n))
print(memory_usage(t, max_usage=True))

"""
Explanation: ```{admonition} Ejercicio
:class: tip
Perfilar la regla de Simpson compuesta que se implementó en el primer ejercicio con line_profiler utilizando su import en una máquina de AWS con las mismas características que la que se presenta en esta nota.
```
Perfilamiento: medición de uso de memoria en Python
Si bien las computadoras de hoy en día tienen una gran cantidad de RAM es importante que las aplicaciones no utilicen la totalidad pues en ese caso se tendrá una penalización en el performance de la aplicación al utilizar virtual memory, ver espacio_de_intercambio.
memory_profiler
Se ejecuta desde la línea de comandos, con un comando de magic o realizando import. Al instalar memory_profiler se incluyen dos comandos de magic: %memit y %mprun. %mprun es similar a line_profiler al analizar línea por línea el uso de memoria.
En el caso de import regresa una lista de valores de uso de memoria en MiB medidas cada cierto interval (argumento de memory_usage). En lo siguiente se pide que se regrese el máximo uso de memoria.
End of explanation
"""


start_mem = memory_usage(max_usage=True)
res = memory_usage(t, max_usage=True, retval=True)
print('start mem', start_mem)
print('max mem', res[0])
print('used mem', res[0]-start_mem)
print('fun output', res[1])

"""
Explanation: Y podemos ejecutar las siguientes líneas de acuerdo a how-do-i-profile-memory-usage-in-python para consistencia con los resultados que se obtendrán con %memit y %mprun:
End of explanation
"""


%load_ext memory_profiler

"""
Explanation: Para el uso de %memit y %mprun se carga la extensión de jupyter:
End of explanation
"""


%memit #how much RAM this process is consuming

%memit -c Rcf(f, a, b, n)

%memit -c Rcf_2(f, a, b, n)

%memit -c Rcf_3(f_np, a, b, n)

"""
Explanation: %memit devuelve el pico de memoria usada en una celda de un jupyter notebook y utiliza las mismas ideas para reportar las mediciones que %timeit:
End of explanation
"""


%%file Rcf_memory_profiler.py

import math

from pytest import approx
from scipy.integrate import quad
from memory_profiler import profile

@profile #decorate the functions you want to profile with memory_profiler 
def Rcf(f,a,b,n):
    """
    Compute numerical approximation using rectangle or mid-point
    method in an interval.
    Nodes are generated via formula: x_i = a+(i+1/2)h_hat for
    i=0,1,...,n-1 and h_hat=(b-a)/n
    Args:
    
        f (float): function expression of integrand.
        
        a (float): left point of interval.
        
        b (float): right point of interval.
        
        n (int): number of subintervals.
        
    Returns:
    
        sum_res (float): numerical approximation to integral
            of f in the interval a,b
    """
    h_hat = (b-a)/n
    nodes = [a+(i+1/2)*h_hat for i in range(n)]
    sum_res = 0
    for node in nodes:
        sum_res = sum_res+f(node)
    return h_hat*sum_res

if __name__ == "__main__":
    a = 0
    b = 1
    f = lambda x: math.exp(-x**2)
    obj, err = quad(f, a, b)
    print("objetivo: {:0.6e}".format(obj))
    n = 10**6
    res = Rcf(f,a,b,n)
    print("aproximación: {:0.6e}".format(res))
    print(res == approx(obj))
    

import math

a = 0
b = 1
f = lambda x: math.exp(-x**2)
n = 10**6

from Rcf_memory_profiler import Rcf

%mprun -f Rcf Rcf(f, a, b, n)

"""
Explanation: ```{margin}
Ver glossary: decorator para la definición de un decorator.
```
```{admonition} Comentario
Si se utiliza %mprun debe crearse script.py y utilizar un decorator.
```
End of explanation
"""


%%bash
python3 Rcf_memory_profiler.py

"""
Explanation: El output de %mprun contiene:


La primera columna hace referencia a la línea de código.


Mem usage es la memoria usada por Python después que tal línea ha sido ejecutada. En el caso anterior se tienen de inicio una cantidad positiva de MiB para cargar el script script.py, realizar imports.


Increment representa la diferencia en memoria de la línea actual respecto a la línea anterior. Las cantidades que se muestran en esta columna no necesariamente indican que en memoria se alojan tales cantidades.


Ocurrences contador del número de veces que se ejecuta tal línea.


Line contents el código a perfilar.


Ejemplo ejecutado desde la línea de comandos:
End of explanation
"""


%%file heapy_ex_1
import math

from guppy import hpy

hp = hpy()
hp.setrelheap() #Everything allocated before this call will not be in the objects you get later.
a = 0
b = 1
n = 10**6
f = lambda x: math.exp(-x**2)
h_hat = (b-a)/n
h = hp.heap()
print(h)

%%bash
python3 heapy_ex_1

"""
Explanation: ```{admonition} Comentarios

De FAQ

Q: How accurate are the results ?
A: This module gets the memory consumption by querying the operating system kernel about the amount of memory the current process has allocated, which might be slightly different from the amount of memory that is actually used by the Python interpreter. Also, because of how the garbage collector works in Python the result might be different between platforms and even between runs.

También memory_profiler contiene un script de nombre mprof que permite monitorear uso de memoria de un proceso (que puede ser un Python script u otro) y también soporta el monitoreo de child processes que fueron spawn de un proceso en un contexto de multiprocessing-contexts-and-start-methods, ver memory_profiler#time-based-memory-usage

```
```{admonition} Ejercicio
:class: tip
Perfilar la regla de Simpson compuesta que se implementó en el primer ejercicio con memory_profiler utilizando su import en una máquina de AWS con las mismas características que la que se presenta en esta nota.
```
heapy
Python utiliza un memory manager para la administración de la memoria dinámica heap. Aunque Python maneja la memoria automáticamente pueden darse situaciones en las que referencias a objetos ya no sean necesarias o bien tener memory leaks. En tal caso podría llamarse al garbage collector (administrado por el memory manager) de forma manual y usar la herramienta de heapy para monitorear la memoria heap.
```{margin}
Ver what-is-a-memory-heap, Python: memory management, does-python-have-a-stack-heap-and-how-is-memory-managed.
```
La herramienta heapy nos ayuda a obtener información acerca de los objetos en la memoria heap. Por ejemplo podemos revisar el número y tamaño de cada objeto que está en la heap. También nos ayuda a encontrar memory leaks.
heapy funciona escribiendo líneas para llamar a la función heap que encierren las secciones de código que queremos monitorear. Subsecuentes llamados a tal función nos proporciona información en una forma incremental.
End of explanation
"""


%%file heapy_ex_1_b
import math

from guppy import hpy

hp = hpy()
hp.setrelheap() #Everything allocated before this call will not be in the objects you get later.
a = 0
b = 1
n = 10**6
f = lambda x: math.exp(-x**2)
h_hat = (b-a)/n
h = hp.heap()
print(h)
print("index 1:")
print(h[1])

%%bash
python3 heapy_ex_1_b

"""
Explanation: El output anterior muestra:


Partition of a set of ... objects que da el número de objetos analizados y clasificados de acuerdo al Kind.


Total size tamaño de la heap por cada llamado a la función heap del monitoreo. Ayuda a las comparaciones del monitoreo subsecuentes.


Index para enumerar cada objeto en la heap. Se puede utilizar para obtener información utilizando índices, por ejemplo h[0] o h[1]:
End of explanation
"""


%%file heapy_ex_1_2
import math

from guppy import hpy

hp = hpy()
hp.setrelheap() #Everything allocated before this call will not be in the objects you get later.
a = 0
b = 1
n = 10**6
f = lambda x: math.exp(-x**2)
h_hat = (b-a)/n
h = hp.heap()
print(h)
print("index 1:")
print(h[1])
nodes = [a+(i+1/2)*h_hat for i in range(n)]
h = hp.heap()
print("After creating nodes")
print(h)    

%%bash
python3 heapy_ex_1_2

%%file heapy_ex_1_3
import math

from guppy import hpy

hp = hpy()
hp.setrelheap() #Everything allocated before this call will not be in the objects you get later.
a = 0
b = 1
n = 10**6
f = lambda x: math.exp(-x**2)
h_hat = (b-a)/n
h = hp.heap()
print(h)
print("index 1:")
print(h[1])
nodes = [a+(i+1/2)*h_hat for i in range(n)]
h = hp.heap()
print("After creating nodes")
print(h)    
sum_res = 0
for node in nodes:
    sum_res = sum_res+f(node)
h = hp.heap()
print("After sum")
print(h)

%%bash
python3 heapy_ex_1_3

"""
Explanation: Count un contador del número de objetos en el reporte de acuerdo al Kind. Enseguida se muestra su porcentaje de acuerdo al número de objetos analizados.


Size el tamaño que ocupa en la heap de acuerdo al Total size. Enseguida se muestra su porcentaje correspondiente.


Cumulative muestra por cada objeto analizado de acuerdo al Kind de lo que acumulan en la heap. El último renglón del reporte equivale al Total size. Enseguida se muestra su porcentaje correspondiente, el último renglón equivale al 100%.
End of explanation
"""


%%file heapy_ex_2
import math

from guppy import hpy

hp = hpy()
hp.setrelheap() #Everything allocated before this call will not be in the objects you get later.
a = 0
b = 1
n = 10**6
f = lambda x: math.exp(-x**2)
h_hat = (b-a)/n
h = hp.heap()
print(h)
print("index 1:")
print(h[1])
nodes = (a+(i+1/2)*h_hat for i in range(n))
h = hp.heap()
print("After creating nodes")
print(h)    
sum_res = 0
for node in nodes:
    sum_res = sum_res+f(node)
h = hp.heap()
print("After sum")
print(h)

%%bash
python3 heapy_ex_2

"""
Explanation: Utilizando generators:
End of explanation
"""


%%file heapy_ex_3
import numpy as np

from guppy import hpy

hp = hpy()
hp.setrelheap() #Everything allocated before this call will not be in the objects you get later.
a = 0
b = 1
n = 10**6
f_np = lambda x: np.exp(-x**2)
h_hat = (b-a)/n
h = hp.heap()
print(h)
print("index 1:")
print(h[1])
aux_vec = np.linspace(a, b, n+1)
nodes= (aux_vec[:-1]+aux_vec[1:])/2
h = hp.heap()
print("After creating nodes")
print(h)    
res = h_hat*np.sum(f_np(nodes))
print("After sum")
print(h)

%%bash
python3 heapy_ex_3

"""
Explanation: Utilizando funcionalidad de NumPy:
End of explanation
"""


c = 0
n = 5
x = [-1]*n
y = [1.5]*n

for i in range(n):
    c += x[i]*y[i]

print(c)

"""
Explanation: ```{admonition} Ejercicio
:class: tip
Perfilar la regla de Simpson compuesta que se implementó en el primer ejercicio con heapy en una máquina de AWS con las mismas características que la que se presenta en esta nota..
```
Basic Linear Algebra Subprograms (BLAS)
De acuerdo a la especificación BLAS, hay diferentes niveles dependiendo del número de datos y cantidad de operaciones a realizar que involucran tales datos.
Operación del producto interno estándar o producto punto (nivel 1 de BLAS)
Consideramos $x,y \in \mathbb{R}^n$. El producto punto entre $x$ y $y$ es $c = x^Ty = \displaystyle \sum_{i=1}^n x_iy_i$. 
Ejemplo y algoritmo
End of explanation
"""


alpha = 2
n = 5
x = [-2]*n
y = [0]*n

for i in range(n):
    y[i] += alpha*x[i]

print(y)

"""
Explanation: ````{admonition} Comentarios


El producto punto de dos $n$-vectores involucran $n$ multiplicaciones y $n$ sumas para un total de $2n$ operaciones o floating point operations per second (flops). Usamos la notación $\mathcal{O}(\cdot)$ para escribir que el producto punto es $\mathcal{O}(n)$ y se lee "de orden $n$ o proporcional a $n$" para indicar que la cantidad de trabajo tiene un comportamiento lineal con la dimensión $n$. También tal cantidad de trabajo opera sobre una cantidad lineal de datos.


Los flops que realiza un algoritmo es una forma de cuantificar el volumen de trabajo asociado con un cálculo. Un flop es una operación de punto flotante: suma, multiplicación o división. Por ejemplo, en la línea:


python
C[i][j] = C[i][j] + A[i][k]*B[k][j]
se realizan $2$ flops. Los flops sólo representan una componente para categorizar a los algoritmos de acuerdo al trabajo que realizan, otras componentes son la transferencia o movimientos de datos, data movement/motion, ejecución secuencial o en paralelo y el data locality y data reuse que realizan.

En LAPACK encontramos sdot, ddot, cdotu, zdotu para descripción de las funciones/subrutinas escritas en Fortran del producto punto en los casos de precisión simple, doble, números complejos precisión simple y doble respectivamente.

````
Operación saxpy (nivel 1 de BLAS)
Consideramos $\alpha \in \mathbb{R}, x,y \in \mathbb{R}^n$. El nombre lo recibe por scalar alpha x plus y. En LAPACK saxpy se escribe en forma update:
$$y=\alpha x + y \therefore y_i = \alpha x_i + y_i \forall i=1,...,n$$
```{admonition} Comentarios


El símbolo "=" no se utiliza como igualdad de expresiones sino un assignment statement.


También encontramos en LAPACK caxpy, zaxpy y daxpy para el caso complejo precisión simple y doble y para números en doble precisión respectivamente.


Ésta operación realiza un trabajo de $\mathcal{O}(n)$ sobre una cantidad de datos $\mathcal{O}(n)$.


```
Ejemplo y algoritmo
End of explanation
"""


alpha = 2
n = 5
x = [-2]*n
y = [3,4,-1,0,1]

for i in range(n):
    y[i] += alpha*x[i]

print(y)

"""
Explanation: o en una forma update:
End of explanation
"""


m = 2
n = 5
A = [[1.2]*n if i%2==0 else [1]*n for i in range(m)]

print(A)

"""
Explanation: ```{admonition} Comentario
La operación de producto punto y saxpy son algoritmos catalogados como de nivel BLAS 1, ver BLAS: Basic Linear Algebra Subprograms. Éstos algoritmos se caracterizan por involucrar una cantidad de trabajo lineal sobre una cantidad lineal de datos. Ver level 1 para más ejemplos de este tipo de algoritmos.
```
Operación de multiplicación matriz-vector (nivel 2 de BLAS)
Consideramos $A \in \mathbb{R}^{m \times n}, x \in \mathbb{R}^n, y \in \mathbb{R}^m$. La operación $y = y + Ax$ es una operación generalizada saxpy, por ello se denomina gaxpy pero en LAPACK podemos encontrarla con nombres como sgemv, dgemv, cgemv y zgemv para los casos de precisión simple, doble, números complejos precisión simple y doble respectivamente. Hay diferentes formas de visualizar y escribir el algoritmo de multiplicación matriz-vector. Por ejemplo para una matriz $A$ con entradas:
End of explanation
"""


x = [2]*n
y = [0]*m
for i in range(m):
    for j in range(n):
        y[i] += A[i][j]*x[j]


print(y)

"""
Explanation: se tiene:
Ejemplo y algoritmo gaxpy row oriented
End of explanation
"""


x = [2]*n
y = [-1]*m
for i in range(m):
    for j in range(n):
        y[i] += A[i][j]*x[j]


print(y)

"""
Explanation: Si $y$ tiene valores distintos de $0$, se realiza un update:
End of explanation
"""


x = [2]*n
y = [0]*m
for j in range(n):
    for i in range(m):
        y[i] += A[i][j]*x[j]

print(y)

"""
Explanation: ````{admonition} Comentarios

En la versión row oriented del algoritmo gaxpy, el inner loop realiza productos punto entre el $i$-ésimo renglón de $A$ y el vector $x$. Se realizan $m$ productos punto $A[i,:]^Tx$

```python
for i in range(m):
    y[i] += A[i,:]*x #producto punto
```
donde: $A[i,:]$ es el $i$-ésimo renglón de $A$. Así podemos reescribir de forma más compacta este algoritmo.
Sin embargo como hemos visto en Python con su implementación más común CPython, no es posible realizar tal indexado pues resulta en el error:
```python
TypeError                                 Traceback (most recent call last)
<ipython-input-14-8ae2a3e9d101> in <module>
      1 for i in range(m):
----> 2     y[i] += A[i,:]*x #producto punto
      3 
TypeError: list indices must be integers or slices, not tuple
```
a menos que incorporemos alguna paquetería que permita la vectorización y el uso de índices para extracción de columnas (o renglones) de $A$ como NumPy.
````
```{admonition} Observación
:class: tip
Obsérvese que el acceso a la matriz $A$ del algoritmo gaxpy row oriented es por renglón, de ahí el nombre.
```
También puede escribirse al algoritmo gaxpy en una forma orientada por columnas.
Ejemplo y algoritmo gaxpy column oriented
Este algoritmo ayuda a visualizar al producto matriz-vector como una combinación lineal de las columnas de $A$:
$$Ax = \displaystyle \sum_{j=1}^n A_jx_j$$
con $A_j$ la $j$-ésima columna de $A$.
End of explanation
"""


x = 2*np.ones(n)
y = np.zeros(m)

print(x)

print(y)

A = np.array([[1.2,1.2,1.2,1.2,1.2],[1,1,1,1,1]])

print(A)

for j in range(n):
    y += A[:,j]*x[j]

print(y)

"""
Explanation: ````{admonition} Comentarios 


El algoritmo de multiplicación matriz-vector (versión row o column oriented) involucra $\mathcal{O}(mn)$ operaciones o una cantidad cuadrática de trabajo, que podemos entender como "si duplicamos cada dimensión de $A$ entonces la cantidad de trabajo se incrementa por un factor de $4$". Tal número de operaciones trabajan sobre una matriz o sobre una cantidad cuadrática de datos. A los algoritmos que realizan una cantidad cuadrática de trabajo sobre una cantidad cuadrática de datos se les cataloga de nivel BLAS 2. Ver level 2 para más ejemplos de algoritmos en el álgebra lineal en esta categoría.


En el algoritmo gaxpy column oriented el acceso a la matriz $A$ es por columna.


La versión column oriented se puede analizar desde el punto de vista puramente algorítmico como un intercambio entre las líneas con los índices $i$ y $j$ de cada loop y un acceso a los datos de la matriz por columna. O bien, se puede analizar desde el álgebra lineal indicando que el vector $y$ está en el espacio generado por las columnas de $A$ y cuyas coordenadas son dadas por las entradas del vector $x$:


<img src="https://dl.dropboxusercontent.com/s/6a2b7rjs4a71sni/combinacion_lineal_columnas_A.png?dl=0" heigth="700" width="700">

Obsérvese que el inner loop de la versión column oriented en gaxpy es un saxpy en la que el escalar está dado por una entrada de $x$. Esto lo podemos escribir de forma explícita definiendo $A[:,j]$ a la $j$-ésima columna de $A$ por lo que $A = [A[:,1] | A[:,2] | \dots | A[:,n]]$, entonces:

python
x = [2]*n
y = [0]*m
for j in range(n):
    y += A[:,j]*x[j]
Sin embargo como hemos visto en Python con su implementación más común CPython, no es posible realizar tal indexado pues resulta en el error:
```python
TypeError                                 Traceback (most recent call last)
<ipython-input-17-93f449af9194> in <module>
      2 y = [0]m
      3 for j in range(n):
----> 4     y += A[:,j]x[j]
TypeError: list indices must be integers or slices, not tuple
```
a menos que incorporemos alguna paquetería que permita la vectorización y el uso de índices para extracción de columnas (o renglones) de $A$ como NumPy.
````
End of explanation
"""


x = 2*np.ones(n)
y = np.zeros(m)
A = np.array([[1.2,1.2,1.2,1.2,1.2],[1,1,1,1,1]])

for i in range(m):
    y[i] += A[i,:].dot(x)

print(y)

"""
Explanation: El algoritmo gaxpy row oriented puede escribirse de forma más compacta haciendo uso de la definición de producto punto estándar: $x^Ty$ para dos vectores columna $x$ y $y$. En el caso de una matriz $A$ se tiene:
python
for i in range(m)
    y[i] += A[i,:]^T*x
donde: $A[i,:]$ es el $i$-ésimo renglón de $A$. En Python:
End of explanation
"""


print(np.show_config())

"""
Explanation: Operación de multiplicación matriz-matriz (nivel 3 de BLAS)
El algoritmo de multiplicación de matrices $C = C + AB$ con $A \in \mathbb{R}^{m \times r}, B \in \mathbb{R}^{r \times n}, C \in \mathbb{R}^{m \times n}$ se cataloga como de nivel 3 de BLAS pues realiza una cantidad de trabajo cúbica sobre una cantidad cuadrática de datos. Ver level 3 para más ejemplos de algoritmos en el álgebra lineal en esta categoría.
En LAPACK encontramos al algoritmo de multiplicación matricial con nombres como sgemm, dgemm y cgemm, zgemm para los casos de precisión simple, doble, números complejos precisión simple y doble respectivamente.
El algoritmo de multiplicación de matrices puede escribirse en diferentes versiones. Por ejemplo la versión  i,j,k es:
python
for i in range(m):
    for j in range(n):
        for k in range(r):
            C[i][j]+=A[i][k]*B[k][j]
y la versión j,k,i es:
python
for j in range(n):
    for k in range(r):
        for i in range(m):
            C[i][j]+=A[i][k]*B[k][j]
```{admonition} Comentarios


Cualquiera de las versiones (hay $3!$ de ellas) involucra una cantidad de trabajo del orden $\mathcal{O}(mnr)$, la cual es cúbica. Es posible interpretar ésta cantidad de trabajo con una frase del tipo "si duplicamos cada dimensión de $A$ entonces la cantidad de trabajo se incrementa por un factor de $8$".


Distintas versiones tienen distinto patrón de acceso a los datos de $A, B$ y $C$. Por ejemplo, para la variante i,j,k el inner loop realiza un producto punto que requiere acceso a renglones de $A$ y columnas de $B$ en una forma izquierda-derecha y arriba-abajo respectivamente. La variante j,k,i involucra una operación saxpy en el inner loop y acceso por columnas a la matriz $C$ y a $A$. Un resúmen de lo anterior se presenta en el siguiente cuadro:


<img src="https://dl.dropboxusercontent.com/s/rw81crmo1b7fjvb/tabla_con_versiones_mult_matrices.png?dl=0" heigth="450" width="450">
```
Otras versiones de la multiplicación de matrices apoyándose de las formas (i,j,k), (j,k,i),...
Si utilizamos la versión i,j,k podemos reescribir la multiplicación de matrices en una forma:
python
for i in range(m):
    for j in range(n)
        C[i][j] += A[i,:].dot(B[:,j])
y hemos vectorizado el inner loop para usar productos punto que es una operación nivel 1 de BLAS.
Además NumPy nos provee de funcionalidad para reescribir ésta forma de productos punto en una como sigue:
python
for i in range(m):
    C[i,:] += A[i,:]@B
```{admonition} Comentario
En esta reescritura se realizan $m$ operaciones gaxpy de nivel 2 en el loop. Esta versión debe ser más rápida que la que utiliza operaciones de nivel 1 de BLAS.
```
(OPERTIPICASALGLINEALYFLOPS)=
Operaciones típicas del álgebra lineal y los flops que realizan
A continuación se muestra una tabla que presenta el número de flops realizados por distintas operaciones del álgebra lineal:
<img src="https://dl.dropboxusercontent.com/s/rvqkokicaqkwrif/tabla_con_flops_para_operaciones_alg_lineal.png?dl=0" heigth="550" width="550">
```{admonition} Comentario
Por la tabla anterior, gaxpy (nivel 2 de BLAS) realiza más flops que un producto punto. Aún así, debe cumplirse que la versión de gaxpy sea más rápida que la versión con productos punto (nivel 1 de BLAS) para la multiplicación de matrices pues aprovecha mejor el data reuse y data locality. Esto se logra pues gaxpy disminuye el número de cache misses y por tanto el tráfico hacia y desde el caché (data movement/motion).
```
Métodos para resolver sistemas de ecuaciones lineales (SEL) o factorizaciones matriciales (FM) y los flops que realizan
Ver {ref}algoritmos para resolver sistemas de ecuaciones lineales &lt;ALGSEL&gt; para ejemplos de métodos directos o basados en FM e iterativos que resuelven SEL.
Los métodos iterativos no realizan FM. Hay características de la matriz (por ejemplo simetría) que no explotan en su implementación pero se benefician de otras (por ejemplo si $A$ es sparse o rala: muchos ceros). Otras características de tales métodos son: uso de menos memoria que los directos o basados en FM, convergen bajo condiciones simples de verificar y son sencillos de implementar.
Los métodos directos o basados en FM en general consumen más memoria que los iterativos, tienen el potencial de utilizar operaciones de nivel BLAS $3$, siempre encuentran solución si la matriz del sistema es no singular y en general requieren mejor dominio de la programación para su implementación que los iterativos. Una vez calculada la FM de la matriz $A$ se utilizan métodos para resolver SEL que resultan de la FM (por ejemplo usar FM y luego sustitución hacia delante y hacia atrás para resolver SEL triangulares).
Nombremos a los métodos o algoritmos numéricos anteriores nuestros métodos básicos del álgebra lineal para resolver SEL o FM (MBAL para SEL o FM). Cada uno puede utilizarse para resolver SEL, FM o para otros propósitos. Por ejemplo la factorización QR también se utiliza en el cálculo de eigenvalores y eigenvectores en el algoritmo QR.
Aquí se presenta una tabla que muestra el número de flops si se desea resolver $Ax=b$ con $A \in \mathbb{R}^{n \times n}$ de algunos MBAL para SEL o FM:
<img src="https://dl.dropboxusercontent.com/s/7i01uxbxyadhuh1/flops_MBAL_ABAL.png?dl=0" heigth="250" width="250">
```{admonition} Comentarios


Es fácil verificar que resolver $Ax = b$ calculando la inversa (por ejemplo np.linalg.inv(A)@b) resultan en aproximadamente $\frac{8n^3}{3}$ flops. 


Los MBAL para SEL o FM sirven de apoyo para resolver SEL por bloques. 
```


Blocking algorithms para multiplicación de matrices
El data reuse y el data locality se incrementa al utilizar operaciones de nivel BLAS mayores. También disminuye el data movement/motion hacia y desde el caché y por tanto tiene menos impacto el bottleneck de Von Neumann en el performance.
Por lo anterior hay algoritmos que se diseñan con el objetivo de aprovechar lo más posible el data reuse y el data locality y se nombran blocking algorithms. Entre los blocking algorithms encontramos a los que trabajan con matrices por bloques pues son más eficientes al utilizar un nivel $3$ de BLAS. Una matriz $A \in \mathbb{R}^{m \times n}$ por bloques se puede escribir como:
<img src="https://dl.dropboxusercontent.com/s/s77weq74yoi2rfc/matrix_A_by_blocks.png?dl=0" heigth="300" width="300">
donde: $m_1 + \dots + m_q = m$, $n_1 + \dots + n_r = n$. Con esta definición se llama a la matriz $A$, una matriz por bloques de tamaño $q \times r$.
```{admonition} Comentarios


Hay que tener en cuenta que trabajar con matrices por bloques utiliza mayor cantidad de memoria que un algoritmo que opera a un nivel escalar esto puede ser benéfico o no dependiendo del problema (por ejemplo simetría de una matriz puede aprovecharse en este caso), máquina o arquitectura en la que se esté trabajando.


Trabajar con blocking algorithms permite simplificar notación matemática. Además, su escritura ayuda a pensar en una implementación con cómputo en paralelo.


```
Algunas operaciones en matrices por bloques
Entre las operaciones que son posibles realizar a un nivel por bloques se encuentran:
Multiplicación por escalares
$$\begin{array}{l}
\mu \left[ \begin{array}{cc}
A_{11} & A_{12}\
A_{21} & A_{22}\
A_{31} & A_{32}
\end{array}
\right] = 
\left[\begin{array}{cc}
\mu A_{11} & \mu A_{12}\
\mu A_{21} & \mu A_{22}\
\mu A_{31} & \mu A_{32}
\end{array}
\right] 
\end{array}
$$
con $\mu \in \mathbb{R}$.
Transposición
$$\begin{array}{l}
\left[ \begin{array}{cc}
A_{11} & A_{12}\
A_{21} & A_{22}\
A_{31} & A_{32}
\end{array}
\right]^T = 
\left[\begin{array}{ccc}
A_{11}^T & A_{21}^T & A_{31}^T\
A_{12}^T & A_{22}^T & A_{32}^T\
\end{array}
\right] 
\end{array}
$$
Suma
$$\begin{array}{l}
\left[ \begin{array}{cc}
A_{11} & A_{12}\
A_{21} & A_{22}\
A_{31} & A_{32}
\end{array}
\right] + 
\left[\begin{array}{cc}
B_{11} & B_{12}\
B_{21} & B_{22}\
B_{31} & B_{32}
\end{array}
\right] =
\left[\begin{array}{cc}
A_{11} + B_{11} & A_{12} + B_{12}\
A_{21} + B_{21} & A_{22} + B_{22}\
A_{31} + B_{31} & A_{32} + B_{32}
\end{array}
\right]
\end{array}
$$
Multiplicación
$$\begin{array}{l}
\left[ \begin{array}{cc}
A_{11} & A_{12}\
A_{21} & A_{22}\
A_{31} & A_{32}
\end{array}
\right] \cdot
\left[\begin{array}{cc}
B_{11} & B_{12}\
B_{21} & B_{22}\
\end{array}
\right] =
\left[\begin{array}{cc}
A_{11}B_{11} + A_{12}B_{21} & A_{11}B_{12} + A_{12}B_{22}\
A_{21}B_{11} + A_{22}B_{21} & A_{21}B_{12} + A_{22}B_{22}\
A_{31}B_{11} + A_{32}B_{21} & A_{31}B_{12} + A_{32}B_{22}
\end{array}
\right]
\end{array}
$$
```{admonition} Comentario
El número de columnas de los bloques $A_{11}, A_{21}, A_{31}$ deben coincidir con el número de renglones de $B_{11}, B_{12}$. Así como con los bloques $A_{12}, A_{22}, A_{32}$ y $B_{21}, B_{22}$.
```
Considérese que se particionan a las matrices $A, B, C$ como sigue:
<img src="https://dl.dropboxusercontent.com/s/j0p89i19vf5r88k/A_B_C_matrices_by_blocks.png?dl=0" heigth="400" width="400">
Esto es, $A,B,C$ son $N \times N$ matrices por bloques con $\ell \times \ell$ bloques. Entonces, para índices $\alpha = 1,2,\dots, N$ y $\beta = 1,2,\dots,N$ se tiene que el bloque $\alpha \beta$ de $C$ se obtiene: $C_{\alpha \beta} = \displaystyle \sum_{\gamma=1}^N A_{\alpha \gamma} B_{\gamma \beta}$. 
El algoritmo una vez definido el bloque de tamaño $\ell$ es:
```python
for alpha in np.arange(N):
    i = np.arange((alpha - 1)l + 1,alphal + 1) #TODO: check if index is well defined
    for beta in np.arange(N):
        j = np.arange((beta - 1)l + 1, betal + 1) #TODO: check if index is well defined
        for gamma in np.arange(N):
            k = np.arange((gamma - 1)l + 1, gammal + 1) #TODO: check if index is well defined
            C[i][:,j]+= A[i][:,k]*B[k][:,j]
```
```{admonition} Ejercicio
:class: tip
Implementar la multiplicación por bloques en Python y R. 
```
Métodos por bloques para resolver sistemas de ecuaciones lineales
Eliminación por bloque
Este método consiste en eliminar un subconjunto de variables y resolver un sistema más pequeño. Si el sistema más pequeño asociado a la submatriz puede resolverse con alguno de los MBAL para SEL o FM anteriores entonces este método puede tener más eficiencia que uno que no trabaja por bloques.
Considérese un sistema de la forma: $Ax=b$ escrito como:
$$\begin{array}{l}
\left[ \begin{array}{cc}
A_{11} & A_{12}\
A_{21} & A_{22}
\end{array}
\right] \cdot 
\left[\begin{array}{c}
x_1\
x_2
\end{array}
\right] =
\left[\begin{array}{c}
b_1\
b_2
\end{array}
\right]
\end{array}
$$
con $x_1 \in \mathbb{R}^{n_1}, x_2 \in \mathbb{R}^{n_2}, A_{11} \in \mathbb{R}^{n_1 \times n_1}, A_{22} \in \mathbb{R}^{n_2 \times n_2}, b_1 \in \mathbb{R}^{n_1}, b_2 \in \mathbb{R}^{n_2}$.
El método de eliminación por bloques consiste en que si $A_{11}$ es invertible, entonces se puede eliminar $x_1$ de las ecuaciones como sigue:


De la primer ecuación por bloques: $A_{11}x_1 + A_{12}x_2 = b_1$ se obtiene: $x_1 = A_{11}^{-1}(b_1-A_{12}x_2)$.


Sustituyendo esta relación en la segunda ecuación por bloques obtenemos la ecuación reducida: $(A_{22}-A_{21}A_{11}^{-1}A_{12})x_2 = b_2 - A_{21}A_{11}^{-1}b_1$.


```{admonition} Comentarios


Las dos ecuaciones anteriores son equivalentes al sistema original $Ax=b$. 


La matriz $A_{22}-A_{21}A_{11}^{-1}A_{12}$ se le nombra complemento de Schur del bloque $A_{11}$ en $A$ (ver Schur complement) esto es: 


$$S = A_{22}-A_{21}A_{11}^{-1}A_{12}.$$

Se tiene la propiedad siguiente: $S$ es no singular si y sólo si $A$ es no singular.

```
Algoritmo: método de eliminación por bloque

Dados $A$ y $A_{11}$ no singulares
1.Calcular $A_{11}^{-1}A_{12}$ y $A_{11}^{-1}b_1$ teniendo cuidado en no calcular la inversa sino un sistema de ecuaciones lineales:

Para realizar la multiplicación $A_{11}^{-1}b_1$ definimos $y=A_{11}^{-1}b_1$ y por tanto $A_{11}y = b_1$ ($A_{11}$ es no singular). Así, resolvemos para $y$ el sistema anterior y habremos calculado $A_{11}^{-1}b_1$. Similarmente definimos $Y=A_{11}^{-1}A_{12}$ con lo que se tiene $A_{11}Y=A_{12}$. Resolvemos para $Y \in \mathbb{R} ^{n_1 \times n_1}$ y habremos calculado $A_{11}^{-1}A_{12}$.

2.Calcular el complemento de Schur del bloque $A_{11}$ en $A$: $S = A_{22}-A_{21}A_{11}^{-1}A_{12}$. Calcular $ \hat{b} = b_2-A_{21}A_{11}^{-1}b_1$.
3.Resolver $Sx_2 = \hat{b}$.
4.Resolver $A_{11}x_1 = b_1-A_{12}x_2$.

```{admonition} Comentario
En los pasos que se requieren resolver sistemas de ecuaciones lineales se utilizan los MBAL para SEL o FM.
```
Factor-Solve method
El método de factor solve consiste en expresar a la matriz del sistema $Ax = b$ como un producto de la forma: $A = A_1A_2 \cdots A_k$. Si $A$ es no singular entonces: $x = A^{-1}b = A_k^{-1}A_{k-1}^{-1} \cdots A_1^{-1}b$.
De esta forma, $x$ puede calcularse de "derecha a izquierda":
Descripción matemática del algoritmo
$$
\begin{eqnarray}
z_1&:=&A_1^{-1}b \nonumber \
z_2&:=&A_2^{-1}z_1 &= A_2^{-1}A_1^{-1}b \nonumber\
\vdots \nonumber\
z_{k-1}&:=&A_{k-1}^{-1}z_{k-2} &= A_{k-1}^{-1} \cdots A_1^{-1}b \nonumber\
x &:=&A_k^{-1}z_{k-1} &= A_k^{-1} \cdot A_{k-1}^{-1} \cdots A_1^{-1}b \nonumber
\end{eqnarray}
$$
```{admonition} Comentarios


El paso que consiste en factorizar a la matriz $A$ es el paso nombrado factor y el usado en resolver el sistema $A_iz_i = z_{i-1}$ es nombrado solve.


Normalmente $A_i$ tiene una estructura aprovechable (por ejemplo triangular superior/inferior) para resolver el sistema $A_iz_i = z_{i-1}$. Por ejemplo en el caso de la factorización $LU$ de $A$ se cumple $PA = LU$ con $P$ matriz de permutación, $L$ triangular inferior con $1$'s en la diagonal y $U$ triangular superior. Entonces:


Paso 1: encontrar factores $P,L,U$ tales que $PA=LU$.


Paso 2: resolver $Ld=Pb$.


Paso 3: resolver $Ux=d$.




```
Método de eliminación por bloques escrito como un factor-solve method:
El método de block elimination puede escribirse como un factor-solve method:
Paso $1$:
$$\begin{array}{l}
\left[ \begin{array}{cc}
A_{11} & A_{12}\
A_{21} & A_{22}
\end{array}
\right]=
\left[\begin{array}{cc}
A_{11} & 0 \
A_{21} & S
\end{array}
\right]
\left[\begin{array}{cc}
I & A_{11}^{-1}A_{12} \
0 & I
\end{array}
\right]
\end{array}
$$
Paso $2$. Resolver por el método de sustitución hacia delante por bloques:
$$\begin{array}{l}
\left[ \begin{array}{cc}
A_{11} & 0\
A_{21} & S
\end{array}
\right] \cdot
\left[\begin{array}{c}
z_1 \
z_2
\end{array}
\right] = 
\left[\begin{array}{c}
b_1\
b_2
\end{array}
\right]
\end{array}
$$
Paso $3$. Resolver por el método de sustitución hacia atrás por bloques:
$$\begin{array}{l}
\left[ \begin{array}{cc}
I & A_{11}^{-1}A_{12}\
0 & I
\end{array}
\right] \cdot
\left[\begin{array}{c}
x_1 \
x_2
\end{array}
\right] = 
\left[\begin{array}{c}
z_1\
z_2
\end{array}
\right]
\end{array}
$$
```{admonition} Comentario
En los pasos que se requieren resolver sistemas de ecuaciones lineales se utilizan los MBAL para SEL o FM.
```
OpenBLAS y NumPy
En esta sección se sugiere para reproducibilidad de resultados utilizar la AMI y la instancia de AWS ambas descritas al inicio.
Resultados utilizando el script_profiling_and_BLAS.sh
Verificamos que NumPy está usando la implementación de la API de BLAS con OpenBLAS que se obtiene con sudo apt-get install -y libopenblas-dev.
End of explanation
"""


np.random.seed(2020)
m=10**4
r=10**4

A=np.random.rand(m,r)

np.random.seed(2021)
r=10**4
n=10**4

B=np.random.rand(r,n)

"""
Explanation: Multiplicación de matrices
Hacemos una prueba de medición de tiempos:
End of explanation
"""


%timeit -n 1 -r 7 A@B

"""
Explanation: ```{admonition} Observación
:class: tip
Mientras se ejecuta la siguiente celda se sugiere en la terminal ejecutar en la línea de comando htop.
```
End of explanation
"""


sessionInfo()

set.seed(2020)
m<-10**4
r<-10**4
A <- matrix(runif(m*r), m, r)

set.seed(2021)
r<-10**4
n<-10**4
B <- matrix(runif(r*n), r, n)

library(microbenchmark)

"""
Explanation: ```{admonition} Comentarios


La implementación de la API de BLAS con OpenBLAS que se instala mediante sudo apt-get install -y libopenblas-dev y NumPy instalado vía pip install numpy --no-binary numpy en el script script_profiling_and_BLAS.sh, es la misma que la integrada en la instalación de NumPy vía pip install numpy en el script script_no_OpenBLAS.sh.


Más cores reducen el tiempo de cómputo, para ver esto probar lo anterior con máquinas con menos cores.


```
OpenBLAS y librerías de cómputo matricial de R
Resultados utilizando el script_profiling_and_BLAS.sh
Verificamos que las librerías de cómputo matricial de R están usando la implementación de la API de BLAS con OpenBLAS que se obtiene con sudo apt-get install -y libopenblas-dev.
End of explanation
"""


mbk<-microbenchmark(
    A%*%B,
    times=7
    )

print(mbk)

"""
Explanation: ```{admonition} Comentario
El paquete de microbenchmark nos ayuda a la medición de estadísticas para pequeñas piezas de código. Se repetirá ésta medición un número definido de veces.
```
```{admonition} Observación
:class: tip
Mientras se ejecuta la siguiente celda se sugiere en la terminal ejecutar en la línea de comando htop.
```
End of explanation
"""


sessionInfo()

set.seed(2020)
m<-10**4
r<-10**4
A <- matrix(runif(m*r), m, r)

set.seed(2021)
r<-10**4
n<-10**4
B <- matrix(runif(r*n), r, n)

library(microbenchmark)

"""
Explanation: Resultados utilizando el script_no_OpenBLAS.sh
Verificamos que las librerías de cómputo matricial de R no están usando la implementación de la API de BLAS con OpenBLAS que se obtiene con sudo apt-get install -y libopenblas-dev.
End of explanation
"""


mbk<-microbenchmark(
    A%*%B,
    times=1 #just one time as it takes too much time
    )

print(mbk) #after ten minutes didnt finish

"""
Explanation: ```{admonition} Observación
:class: tip
Mientras se ejecuta la siguiente celda se sugiere en la terminal ejecutar en la línea de comando htop.
```
End of explanation
"""


%%file Rcf.py
import math

from pytest import approx
from scipy.integrate import quad
def Rcf(f,a,b,n):
    """
    Compute numerical approximation using rectangle or mid-point
    method in an interval.
    Nodes are generated via formula: x_i = a+(i+1/2)h_hat for
    i=0,1,...,n-1 and h_hat=(b-a)/n
    Args:
    
        f (float): function expression of integrand.
        
        a (float): left point of interval.
        
        b (float): right point of interval.
        
        n (int): number of subintervals.
        
    Returns:
    
        sum_res (float): numerical approximation to integral
            of f in the interval a,b
    """
    h_hat = (b-a)/n
    nodes = [a+(i+1/2)*h_hat for i in range(n)]
    sum_res = 0
    for node in nodes:
        sum_res = sum_res+f(node)
    return h_hat*sum_res
if __name__=="__main__": #añadimos este bloque para ejecución de la función Rcf
    a = 0
    b = 1
    f = lambda x: math.exp(-x**2)
    obj, err = quad(f, a, b)
    print("objetivo: {:0.6e}".format(obj))
    n = 10**6
    res = Rcf(f,a,b,n)
    print("aproximación: {:0.6e}".format(res))
    print(res == approx(obj))

%%bash
/usr/bin/time -p python3 Rcf.py

%%bash
/usr/bin/time -v python3 Rcf.py

"""
Explanation: ````{admonition} Comentarios


La implementación de la API de BLAS con OpenBLAS que se instala mediante sudo apt-get install -y libopenblas-dev y las librerías de cómputo matricial de R que utilizan tal instalación resulta para el algoritmo de multiplicación de matrices ser mucho más rápido que con la implementación de la API de BLAS integrada en la instalación de R.


Más cores reducen el tiempo de cómputo, para ver esto probar lo anterior con máquinas con menos cores.


En R se puede cambiar de implementación de la API de BLAS con el comando: 


sudo update-alternatives --config libblas.so.3-x86_64-linux-gnu
````
Perfilamiento: medición de tiempo en un sistema basado en Unix. Aquí se usa el script_profiling_and_BLAS.sh
/usr/bin/time
Ver también: Wikipedia: Time_Unix
End of explanation
"""


%%file norm_square.py
n = 10**7
vector = list(range(n))
norm = 0
for v in vector:
    norm += v*v

"""
Explanation: ```{admonition} Comentario
Una breve explicación del output se puede encontrar /usr/bin/time. Por ejemplo, para el caso de Major (requiring I/O) nos interesa que sea $0$ pues indica que el sistema operativo tiene que cargar páginas de datos del disco pues tales datos ya no residen en RAM (por alguna razón).
```
```{admonition} Observación
:class: tip
La medición real que regresa /usr/bin/time es el equivalente al Wall time de %time.
```
Perfilamiento: medición de uso de CPU en el sistema operativo GNU/Linux
perf
Ejemplo de cálculo de la norma $2$ al cuadrado de un vector
End of explanation
"""


%%bash
echo "-1" |sudo tee -a /proc/sys/kernel/perf_event_paranoid

%%bash
perf stat -S -e cycles,instructions,cache-references,cache-misses -r 20 python3 norm_square.py

"""
Explanation: ```{margin}
Esta línea de echo es importante para permitir que user ubuntu ejecute perf sin sudo.
```
End of explanation
"""


%%bash
perf stat -S --all-cpus -A -e cycles,instructions,cache-references,cache-misses -r 20 python3 norm_square.py

"""
Explanation: ```{admonition} Comentarios


Con perf se repiten las mediciones utilizando la flag -r. Con la flag -e se enlistan las métricas a calcular. -S llama a sync antes de iniciar la ejecución del programa. -A y --all-cpus para no agregar los conteos para los cores monitoreados. Ver perf stat --help.


En el ejemplo anterior además de los ciclos e instrucciones se calculan los cache references y los cache misses.  Esencialmente el sistema de memoria buscó el número que aparece en cache-references por datos o instrucciones en el caché y de éstas búsquedas, en cache-misses, se repora el número (y porcentaje) de búsquedas fallidas que no estuvieron en memoria. 


Recuérdese que la fragmentación de datos incrementa el número de cache misses.


Se reporta el IPC al lado de instructions con nombre insn per cycle. Un alto IPC indica una alta transferencia de instrucciones de la unidad de memoria hacia la unidad de cómputo y un menor IPC indica más stall cycles, ver pipeline stall. Se debe evaluar si a high rate of instructions indica un high rate of actual work completed (p.ej. un loop tiene una alta tasa de IPC pero no siempre se realiza trabajo útil). Ver sección CPU statistics en la liga. perf también tiene una métrica para medir el número de ciclos que se utilizaron para esperar a ejecutar instrucciones (los cores están stalled). Las métricas son: stalled-cycles-frontend y stalled-cycles-backend.


```
También podemos obtener las estadísticas por core. Obsérvese que la salida anterior de perf calculó el IPC a partir del número de instrucciones y ciclos ahí reportados. El mismo cálculo se realiza a continuación:
End of explanation
"""


%%bash
perf stat -S -r 20 python3 norm_square.py

"""
Explanation: Otras métricas pueden ser obtenidas si ejecutamos perf sólo con la flag -r:
End of explanation
"""


%%bash
perf stat -S --all-cpus -A -r 20 python3 norm_square.py

"""
Explanation: ```{admonition} Comentarios


La métrica de task-clock indica cuántos clock cycles tomó nuestra tarea y se reporta en milisegundos. Obsérvese que también se indica cuántos CPU's fueron utilizados. En el output anterior puede obtenerse no exactamente igual a 1 para CPUs utilized pues el programa no sólo involucra trabajo de CPU sino también de alojamiento de memoria. 


context-switches y cpu-migrations indican cuánto se detuvo el programa para realizar:

operaciones relacionadas con el kernel del sistema (por ejemplo I/O).
ejecución de otras aplicaciones.
alojamiento de la ejecución en un core distinto (lo que ocasiona que haya movimiento de datos o instrucciones hacia el caché de otro core).



la idea es que nuestro programa tenga un número pequeño de éstas métricas.


page-faults se relaciona con el alojamiento de memoria a partir del kernel del sistema operativo. Ocasiona que se detenga la ejecución del programa y por tanto deseamos que ésta métrica no sea muy grande. Ver Page fault.


Para el branching interesa que el número de branch misses sea pequeño.


```
Podemos generar la salida anterior por core:
End of explanation
"""


%%file norm_square_numpy.py
import numpy as np
n = 10**7
vector = np.arange(n)
vector.dot(vector)

%%bash
perf stat -S -e cycles,instructions,cache-references,cache-misses -r 20 python3 norm_square_numpy.py

"""
Explanation: Para contrastar las salidas anteriores con un mejor uso del caché y de los cores se utiliza a continuación el paquete de numpy:
End of explanation
"""


%%bash
perf stat -S --all-cpus -A -e cycles,instructions,cache-references,cache-misses -r 20 python3 norm_square_numpy.py

"""
Explanation: ```{admonition} Comentarios


Obsérvese que es más rápido este programa para el número $n$ que se está utilizando en comparación con el programa anterior sin vectorizar.


El número de instrucciones y el número de cache-misses son menores al reportado con el programa anterior sin vectorizar.


```
Métricas por core:
End of explanation
"""


%%bash
perf stat -S -r 20 python3 norm_square_numpy.py

%%bash
perf stat -S --all-cpus -A -r 20 python3 norm_square_numpy.py

"""
Explanation: ```{admonition} Comentarios


Se observa que la métrica IPC está más uniforme a través de los cores lo que indica un mejor load balancing.


Además, el $\%$ de cache-misses es menor que en el caso no vectorizado.


```
También podemos obtener las salidas anteriores de perf al usar numpy:
End of explanation
"""


np.random.seed(2020)
m = 10**4
n = 10**4
A = np.random.rand(m,n)
file = 'A.txt'
np.savetxt(file,A)

"""
Explanation: perf OpenBLAS: gaxpy row, column oriented (nivel 2)
Utilizaremos perf en lo que sigue para evaluar las métricas de IPC, cache-references, cache-misses y medir el tiempo de ejecución. Además, esto permitirá comparar los tiempos de ejecución del algoritmo gaxpy (nivel BLAS 2), con vectorización y sin vectorización.
End of explanation
"""


%%file mult_matrix_vector.py
m = 10**4
n = 10**4
x = [2.5]*n
y = [0]*m
A = []
file = 'A.txt'
with open(file,'r') as f:
    for l in f:
        A.append([float(k) for k in l.replace('\n','').replace(' ',',').split(',')])      
for j in range(n):
    for i in range(m):
        y[i]+=A[i][j]*x[j]

"""
Explanation: Utilizando listas para construir a la matriz y a los vectores:
End of explanation
"""


%%bash
perf stat -S -e cycles,instructions,cache-references,cache-misses -r 7 python3 mult_matrix_vector.py

"""
Explanation: gaxpy column oriented caso no vectorizado
End of explanation
"""


%%bash
perf stat -S --all-cpus -A -e cycles,instructions,cache-references,cache-misses -r 7 python3 mult_matrix_vector.py

"""
Explanation: Métricas por core:
End of explanation
"""


%%file mult_matrix_vector_numpy.py
import numpy as np
m = 10**4
n = 10**4
x = 2.5*np.ones(n)
y = np.zeros(m)
file = 'A.txt'
A = np.loadtxt(file)
for j in np.arange(n):
    y += A[:,j]*x[j]

%%bash
perf stat -S -e cycles,instructions,cache-references,cache-misses -r 7 python3 mult_matrix_vector_numpy.py

"""
Explanation: gaxpy column oriented caso vectorizando
Vectorizando con NumPy:
End of explanation
"""


%%bash
perf stat -S --all-cpus -A -e cycles,instructions,cache-references,cache-misses -r 7 python3 mult_matrix_vector_numpy.py

"""
Explanation: Métricas por core:
End of explanation
"""


%%file mult_matrix_vector_numpy_row_oriented.py
import numpy as np
m = 10**4
n = 10**4
x = 2.5*np.ones(n)
y = np.zeros(m)
file = 'A.txt'
A = np.loadtxt(file)
for i in np.arange(m):
    y[i] += A[i,:].dot(x)

%%bash
perf stat -S -e cycles,instructions,cache-references,cache-misses -r 7 python3 mult_matrix_vector_numpy_row_oriented.py

%%bash
perf stat -S --all-cpus -A -e cycles,instructions,cache-references,cache-misses -r 7 python3 mult_matrix_vector_numpy_row_oriented.py

"""
Explanation: gaxpy row oriented caso vectorizando
End of explanation
"""


%%file mult_matrix_matrix_numpy_openblas.py
import numpy as np
m = 10**4
r = 10**4
n = 10**4

np.random.seed(2020)
m = 10**4
r = 10**4
A = np.random.rand(m,r)

np.random.seed(2021)
r = 10**4
n = 10**4
B = np.random.rand(r,n)

A@B

%%bash
perf stat -S -e cycles,instructions,cache-references,cache-misses -r 7 python3 mult_matrix_matrix_numpy_openblas.py

%%bash
perf stat -S --all-cpus -A -e cycles,instructions,cache-references,cache-misses -r 7 python3 mult_matrix_matrix_numpy_openblas.py

"""
Explanation: ```{admonition} Comentarios


No se observa una gran diferencia en tiempo de ejecución entre el caso vectorizado y no vectorizado en el nivel 2 de BLAS para la operación gaxpy salvo que en el caso vectorizando se tiene un $\%$ mayor de cache-misses. 


Otra opción que se tiene es utilizar la siguiente implementación ya que NumPy soporta operaciones del tipo matriz-vector sin tener que utilizar un for loop.


```
perf OpenBLAS: nivel 3 multiplicación de matrices
End of explanation
"""


np.random.seed(2020)
m = 10**3
r = 10**3

A = np.random.rand(m,r)
fileA_10_3 = 'A_10_3.txt'
np.savetxt(fileA_10_3,A)

np.random.seed(2021)
m = 10**3
r = 10**3

B = np.random.rand(m,r)
fileB_10_3 = 'B_10_3.txt'
np.savetxt(fileB_10_3,B)

%%file mult_matrix_matrix_numpy_dot_product.py
import numpy as np
m = 10**3
n = 10**3

fileA_10_3 = 'A_10_3.txt'
fileB_10_3 = 'B_10_3.txt'
A = np.loadtxt(fileA_10_3)
B = np.loadtxt(fileB_10_3)
C = np.zeros((m,n))
for i in np.arange(m):
        for j in np.arange(n):
                C[i][j] += A[i,:].dot(B[:,j])

%%bash
perf stat -S --all-cpus -e cycles,instructions,cache-references,cache-misses -r 7 python3 mult_matrix_matrix_numpy_dot_product.py

%%bash
perf stat -S --all-cpus -A -e cycles,instructions,cache-references,cache-misses -r 7 python3 mult_matrix_matrix_numpy_dot_product.py

%%file mult_matrix_matrix_numpy_dot_product_gaxpy.py
import numpy as np
m = 10**3
n = 10**3
fileA_10_3 = 'A_10_3.txt'
fileB_10_3 = 'B_10_3.txt'
A = np.loadtxt(fileA_10_3)
B = np.loadtxt(fileB_10_3)
C = np.zeros((m,n))
for i in np.arange(m):
    C[i,:] = A[i,:]@B

%%bash
perf stat -S --all-cpus -e cycles,instructions,cache-references,cache-misses -r 5 python3 mult_matrix_matrix_numpy_dot_product_gaxpy.py

%%bash
perf stat -S --all-cpus -A -e cycles,instructions,cache-references,cache-misses -r 5 python3 mult_matrix_matrix_numpy_dot_product_gaxpy.py

"""
Explanation: Comparación de versiones de multiplicación de matrices utilizando nivel 1 vs nivel 2 de BLAS
Utilizamos matrices $A, B \in \mathbb{R}^{10^3 \times 10^3}$ pseudoaleatorias:
End of explanation
"""


