.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/plot_benchmark_dot_mul.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_plot_benchmark_dot_mul.py>`
        to download the full example code

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_plot_benchmark_dot_mul.py:


.. _l-example-mul:

Compares matrix multiplication implementations
==============================================

:epkg:`numpy` has a very fast implementation of
matrix multiplication. There are many ways to be slower.

Compared implementations:

* :func:`dmul_cython_omp <td3a_cpp.tutorial.mul_cython_omp.dmul_cython_omp>`
  `code <https://github.com/sdpython/td3a_cpp/blob/master/
  td3a_cpp/tutorial/mul_cython_omp.pyx#L171>`_

.. contents::
    :local:

.. GENERATED FROM PYTHON SOURCE LINES 20-32

.. code-block:: default


    import pprint
    import numpy
    from numpy.testing import assert_almost_equal
    import matplotlib.pyplot as plt
    from pandas import DataFrame, concat
    from td3a_cpp.tutorial.mul_cython_omp import dmul_cython_omp
    from td3a_cpp.tools import measure_time_dim

    dfs = []
    sets = list(range(2, 145, 20))








.. GENERATED FROM PYTHON SOURCE LINES 33-36

numpy mul
+++++++++


.. GENERATED FROM PYTHON SOURCE LINES 36-49

.. code-block:: default


    ctxs = [dict(va=numpy.random.randn(n, n).astype(numpy.float64),
                 vb=numpy.random.randn(n, n).astype(numpy.float64),
                 mul=lambda x, y: x @ y,
                 x_name=n)
            for n in sets]

    res = list(measure_time_dim('mul(va, vb)', ctxs, verbose=1))
    dfs.append(DataFrame(res))
    dfs[-1]['fct'] = 'numpy'
    pprint.pprint(dfs[-1].tail(n=2))






.. rst-class:: sphx-glr-script-out

 .. code-block:: none


      0%|          | 0/8 [00:00<?, ?it/s]
     50%|#####     | 4/8 [00:00<00:00, 16.40it/s]
     75%|#######5  | 6/8 [00:00<00:00, 13.80it/s]
    100%|##########| 8/8 [00:00<00:00,  8.24it/s]
    100%|##########| 8/8 [00:00<00:00,  9.65it/s]
        average  deviation  min_exec  max_exec  ...  number  context_size  x_name    fct
    6  0.000330   0.000001  0.000329  0.000333  ...      50           232     122  numpy
    7  0.000478   0.000002  0.000477  0.000484  ...      50           232     142  numpy

    [2 rows x 9 columns]




.. GENERATED FROM PYTHON SOURCE LINES 50-53

Simple multiplication
+++++++++++++++++++++


.. GENERATED FROM PYTHON SOURCE LINES 53-64

.. code-block:: default


    ctxs = [dict(va=numpy.random.randn(n, n).astype(numpy.float64),
                 vb=numpy.random.randn(n, n).astype(numpy.float64),
                 mul=dmul_cython_omp,
                 x_name=n)
            for n in sets]

    res = list(measure_time_dim('mul(va, vb)', ctxs, verbose=1))
    pprint.pprint(res[-1])






.. rst-class:: sphx-glr-script-out

 .. code-block:: none


      0%|          | 0/8 [00:00<?, ?it/s]
     38%|###7      | 3/8 [00:00<00:00, 19.71it/s]
     62%|######2   | 5/8 [00:01<00:00,  3.33it/s]
     75%|#######5  | 6/8 [00:02<00:01,  1.66it/s]
     88%|########7 | 7/8 [00:05<00:01,  1.13s/it]
    100%|##########| 8/8 [00:09<00:00,  1.95s/it]
    100%|##########| 8/8 [00:09<00:00,  1.17s/it]
    {'average': 0.008107981523498892,
     'context_size': 232,
     'deviation': 2.1798378263439107e-06,
     'max_exec': 0.008112866496667266,
     'min_exec': 0.008105753352865577,
     'number': 50,
     'repeat': 10,
     'x_name': 142}




.. GENERATED FROM PYTHON SOURCE LINES 65-71

Other scenarios
+++++++++++++++

3 differents algorithms, each of them parallelized.
See :func:`dmul_cython_omp
<td3a_cpp.tutorial.mul_cython_omp.dmul_cython_omp>`.

.. GENERATED FROM PYTHON SOURCE LINES 71-87

.. code-block:: default


    for algo in range(0, 2):
        for parallel in (0, 1):
            print("algo=%d parallel=%d" % (algo, parallel))
            ctxs = [dict(va=numpy.random.randn(n, n).astype(numpy.float64),
                         vb=numpy.random.randn(n, n).astype(numpy.float64),
                         mul=lambda x, y: dmul_cython_omp(
                x, y, algo=algo, parallel=parallel),
                x_name=n)
                for n in sets]

            res = list(measure_time_dim('mul(va, vb)', ctxs, verbose=1))
            dfs.append(DataFrame(res))
            dfs[-1]['fct'] = 'a=%d-p=%d' % (algo, parallel)
            pprint.pprint(dfs[-1].tail(n=2))





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    algo=0 parallel=0

      0%|          | 0/8 [00:00<?, ?it/s]
     38%|###7      | 3/8 [00:00<00:00, 18.89it/s]
     62%|######2   | 5/8 [00:01<00:00,  3.30it/s]
     75%|#######5  | 6/8 [00:02<00:01,  1.66it/s]
     88%|########7 | 7/8 [00:05<00:01,  1.12s/it]
    100%|##########| 8/8 [00:09<00:00,  1.93s/it]
    100%|##########| 8/8 [00:09<00:00,  1.17s/it]
        average  deviation  min_exec  ...  context_size  x_name      fct
    6  0.005033   0.000002  0.005031  ...           232     122  a=0-p=0
    7  0.008089   0.000003  0.008085  ...           232     142  a=0-p=0

    [2 rows x 9 columns]
    algo=0 parallel=1

      0%|          | 0/8 [00:00<?, ?it/s]
     50%|#####     | 4/8 [00:00<00:00, 21.80it/s]
     88%|########7 | 7/8 [00:01<00:00,  4.26it/s]
    100%|##########| 8/8 [00:02<00:00,  3.41it/s]
        average  deviation  min_exec  ...  context_size  x_name      fct
    6  0.001276   0.000086  0.001243  ...           232     122  a=0-p=1
    7  0.001861   0.000002  0.001859  ...           232     142  a=0-p=1

    [2 rows x 9 columns]
    algo=1 parallel=0

      0%|          | 0/8 [00:00<?, ?it/s]
     38%|###7      | 3/8 [00:00<00:00, 17.53it/s]
     62%|######2   | 5/8 [00:01<00:00,  3.04it/s]
     75%|#######5  | 6/8 [00:02<00:01,  1.53it/s]
     88%|########7 | 7/8 [00:05<00:01,  1.22s/it]
    100%|##########| 8/8 [00:10<00:00,  2.11s/it]
    100%|##########| 8/8 [00:10<00:00,  1.27s/it]
        average  deviation  min_exec  ...  context_size  x_name      fct
    6  0.005519   0.000001  0.005517  ...           232     122  a=1-p=0
    7  0.008855   0.000003  0.008853  ...           232     142  a=1-p=0

    [2 rows x 9 columns]
    algo=1 parallel=1

      0%|          | 0/8 [00:00<?, ?it/s]
     50%|#####     | 4/8 [00:00<00:00, 20.76it/s]
     88%|########7 | 7/8 [00:01<00:00,  4.36it/s]
    100%|##########| 8/8 [00:02<00:00,  3.44it/s]
        average  deviation  min_exec  ...  context_size  x_name      fct
    6  0.001246   0.000008   0.00124  ...           232     122  a=1-p=1
    7  0.001878   0.000005   0.00187  ...           232     142  a=1-p=1

    [2 rows x 9 columns]




.. GENERATED FROM PYTHON SOURCE LINES 88-93

One left issue
++++++++++++++

Will you find it in :func:`dmul_cython_omp
<td3a_cpp.tutorial.mul_cython_omp.dmul_cython_omp>`.

.. GENERATED FROM PYTHON SOURCE LINES 93-109

.. code-block:: default



    va = numpy.random.randn(3, 4).astype(numpy.float64)
    vb = numpy.random.randn(4, 5).astype(numpy.float64)
    numpy_mul = va @ vb

    try:
        for a in range(0, 50):
            wrong_mul = dmul_cython_omp(va, vb, algo=2, parallel=1)
            assert_almost_equal(numpy_mul, wrong_mul)
            print("Iteration %d is Ok" % a)
        print("All iterations are unexpectedly Ok. Don't push your luck.")
    except AssertionError as e:
        print(e)






.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    Arrays are not almost equal to 7 decimals

    Mismatched elements: 3 / 15 (20%)
    Max absolute difference: 0.98678316
    Max relative difference: 2.10085769e-16
     x: array([[-2.1521276, -1.4024366, -0.8470817,  1.8995522, -0.9867832],
           [-0.1073382,  1.9952894, -0.1321154,  1.6005042,  0.0930345],
           [-1.10652  , -2.1936239, -0.5114178,  0.0623881, -0.3139902]])
     y: array([[-2.1521276, -1.4024366, -0.8470817,  1.8995522,  0.       ],
           [-0.1073382,  1.9952894, -0.1321154,  1.6005042,  0.       ],
           [-1.10652  , -2.1936239, -0.5114178,  0.0623881,  0.       ]])




.. GENERATED FROM PYTHON SOURCE LINES 110-115

Other scenarios but transposed
++++++++++++++++++++++++++++++

Same differents algorithms but the second matrix
is transposed first: ``b_trans=1``.

.. GENERATED FROM PYTHON SOURCE LINES 115-133

.. code-block:: default



    for algo in range(0, 2):
        for parallel in (0, 1):
            print("algo=%d parallel=%d transposed" % (algo, parallel))
            ctxs = [dict(va=numpy.random.randn(n, n).astype(numpy.float64),
                         vb=numpy.random.randn(n, n).astype(numpy.float64),
                         mul=lambda x, y: dmul_cython_omp(
                x, y, algo=algo, parallel=parallel, b_trans=1),
                x_name=n)
                for n in sets]

            res = list(measure_time_dim('mul(va, vb)', ctxs, verbose=2))
            dfs.append(DataFrame(res))
            dfs[-1]['fct'] = 'a=%d-p=%d-T' % (algo, parallel)
            pprint.pprint(dfs[-1].tail(n=2))






.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    algo=0 parallel=0 transposed

      0%|          | 0/8 [00:00<?, ?it/s]
     38%|###7      | 3/8 [00:00<00:00, 25.04it/s]
     75%|#######5  | 6/8 [00:01<00:00,  2.69it/s]
    100%|##########| 8/8 [00:06<00:00,  1.02s/it]
    100%|##########| 8/8 [00:06<00:00,  1.28it/s]
        average  deviation  min_exec  ...  context_size  x_name        fct
    6  0.003414   0.000003  0.003411  ...           232     122  a=0-p=0-T
    7  0.005236   0.000003  0.005234  ...           232     142  a=0-p=0-T

    [2 rows x 9 columns]
    algo=0 parallel=1 transposed

      0%|          | 0/8 [00:00<?, ?it/s]
     50%|#####     | 4/8 [00:00<00:00, 33.63it/s]
    100%|##########| 8/8 [00:01<00:00,  7.00it/s]
    100%|##########| 8/8 [00:01<00:00,  7.94it/s]
        average  deviation  min_exec  ...  context_size  x_name        fct
    6  0.000528   0.000058  0.000504  ...           232     122  a=0-p=1-T
    7  0.000737   0.000003  0.000735  ...           232     142  a=0-p=1-T

    [2 rows x 9 columns]
    algo=1 parallel=0 transposed

      0%|          | 0/8 [00:00<?, ?it/s]
     38%|###7      | 3/8 [00:00<00:00, 25.02it/s]
     75%|#######5  | 6/8 [00:01<00:00,  2.74it/s]
    100%|##########| 8/8 [00:06<00:00,  1.01it/s]
    100%|##########| 8/8 [00:06<00:00,  1.30it/s]
        average  deviation  min_exec  ...  context_size  x_name        fct
    6  0.003318   0.000001  0.003316  ...           232     122  a=1-p=0-T
    7  0.005143   0.000001  0.005142  ...           232     142  a=1-p=0-T

    [2 rows x 9 columns]
    algo=1 parallel=1 transposed

      0%|          | 0/8 [00:00<?, ?it/s]
     50%|#####     | 4/8 [00:00<00:00, 33.17it/s]
    100%|##########| 8/8 [00:00<00:00,  7.15it/s]
    100%|##########| 8/8 [00:00<00:00,  8.10it/s]
        average  deviation  min_exec  ...  context_size  x_name        fct
    6  0.000502   0.000002  0.000501  ...           232     122  a=1-p=1-T
    7  0.000726   0.000003  0.000723  ...           232     142  a=1-p=1-T

    [2 rows x 9 columns]




.. GENERATED FROM PYTHON SOURCE LINES 134-136

Let's display the results
+++++++++++++++++++++++++

.. GENERATED FROM PYTHON SOURCE LINES 136-158

.. code-block:: default


    cc = concat(dfs)
    cc['N'] = cc['x_name']

    fig, ax = plt.subplots(3, 2, figsize=(10, 8), sharex=True, sharey=True)
    ccnp = cc.fct == 'numpy'
    cct = cc.fct.str.contains('-T')
    cca0 = cc.fct.str.contains('a=0')
    cc[ccnp | (~cct & cca0)].pivot(
        'N', 'fct', 'average').plot(logy=True, logx=True, ax=ax[0, 0])
    cc[ccnp | (~cct & ~cca0)].pivot(
        'N', 'fct', 'average').plot(logy=True, logx=True, ax=ax[0, 1])
    cc[ccnp | (cct & cca0)].pivot(
        'N', 'fct', 'average').plot(logy=True, logx=True, ax=ax[1, 0])
    cc[ccnp | (~cct & ~cca0)].pivot(
        'N', 'fct', 'average').plot(logy=True, logx=True, ax=ax[1, 1])
    cc[ccnp | cca0].pivot(index='N', columns='fct', values='average').plot(
        logy=True, logx=True, ax=ax[2, 0])
    cc[ccnp | ~cca0].pivot(index='N', columns='fct', values='average').plot(
        logy=True, logx=True, ax=ax[2, 1])
    fig.suptitle("Comparison of matrix multiplication implementations")



.. rst-class:: sphx-glr-script-out

.. code-block:: pytb

    Traceback (most recent call last):
      File "somewheretd3a_cpp_39_std/td3a_cpp/examples/plot_benchmark_dot_mul.py", line 144, in <module>
        cc[ccnp | (~cct & cca0)].pivot(
    TypeError: pivot() takes 1 positional argument but 4 were given




.. GENERATED FROM PYTHON SOURCE LINES 159-162

The results depends on the machine, its
number of cores, the compilation settings
of :epkg:`numpy` or this module.

.. GENERATED FROM PYTHON SOURCE LINES 162-164

.. code-block:: default


    plt.show()


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** ( 0 minutes  51.066 seconds)


.. _sphx_glr_download_auto_examples_plot_benchmark_dot_mul.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example




    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_benchmark_dot_mul.py <plot_benchmark_dot_mul.py>`

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_benchmark_dot_mul.ipynb <plot_benchmark_dot_mul.ipynb>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_