:orphan:

=================
pentapy Tutorials
=================

In the following you will find Tutorials on how to use pentapy.


Introduction: Solving a pentadiagonal system
============================================

Pentadiagonal systems arise in many areas of science and engineering,
for example in solving differential equations with a finite difference sceme.


Theoretical Background
----------------------

A pentadiagonal system is given by the equation: :math:`M\cdot X = Y`, where
:math:`M` is a quadratic :math:`n\times n` matrix given by:

.. math::

    M=\left(\begin{matrix}d_{1} & d_{1}^{\left(1\right)} & d_{1}^{\left(2\right)} & 0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0\\
    d_{2}^{\left(-1\right)} & d_{2} & d_{2}^{\left(1\right)} & d_{2}^{\left(2\right)} & 0 & \cdots & \cdots & \cdots & \cdots & 0\\
    d_{3}^{\left(-2\right)} & d_{3}^{\left(-1\right)} & d_{3} & d_{3}^{\left(1\right)} & d_{3}^{\left(2\right)} & 0 & \cdots & \cdots & \cdots & 0\\
    0 & d_{4}^{\left(-2\right)} & d_{4}^{\left(-1\right)} & d_{4} & d_{4}^{\left(1\right)} & d_{4}^{\left(2\right)} & 0 & \cdots & \cdots & 0\\
    \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots\\
    \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots\\
    \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots\\
    0 & \cdots & \cdots & \cdots & 0 & d_{n-2}^{\left(-2\right)} & d_{n-2}^{\left(-1\right)} & d_{n-2} & d_{n-2}^{\left(1\right)} & d_{n-2}^{\left(2\right)}\\
    0 & \cdots & \cdots & \cdots & \cdots & 0 & d_{n-1}^{\left(-2\right)} & d_{n-1}^{\left(-1\right)} & d_{n-1} & d_{n-1}^{\left(1\right)}\\
    0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 & d_{n}^{\left(-2\right)} & d_{n}^{\left(-1\right)} & d_{n}
    \end{matrix}\right)

The aim is now, to solve this equation for :math:`X`.


Memory efficient storage
------------------------

To store a pentadiagonal matrix memory efficient, there are two options:

1. row-wise storage:

.. math::

    M_{\mathrm{row}}=\left(\begin{matrix}d_{1}^{\left(2\right)} & d_{2}^{\left(2\right)} & d_{3}^{\left(2\right)} & \cdots & d_{n-2}^{\left(2\right)} & 0 & 0\\
    d_{1}^{\left(1\right)} & d_{2}^{\left(1\right)} & d_{3}^{\left(1\right)} & \cdots & d_{n-2}^{\left(1\right)} & d_{n-1}^{\left(1\right)} & 0\\
    d_{1} & d_{2} & d_{3} & \cdots & d_{n-2} & d_{n-1} & d_{n}\\
    0 & d_{2}^{\left(-1\right)} & d_{3}^{\left(-1\right)} & \cdots & d_{n-2}^{\left(-1\right)} & d_{n-1}^{\left(-1\right)} & d_{n}^{\left(-1\right)}\\
    0 & 0 & d_{3}^{\left(-2\right)} & \cdots & d_{n-2}^{\left(-2\right)} & d_{n-1}^{\left(-2\right)} & d_{n}^{\left(-2\right)}
    \end{matrix}\right)

Here we see, that the numbering in the above given matrix was aiming at the
row-wise storage. That means, the indices were taken from the row-indices of
the entries.


2. column-wise storage:

.. math::

    M_{\mathrm{col}}=\left(\begin{matrix}0 & 0 & d_{1}^{\left(2\right)} & \cdots & d_{n-4}^{\left(2\right)} & d_{n-3}^{\left(2\right)} & d_{n-2}^{\left(2\right)}\\
    0 & d_{1}^{\left(1\right)} & d_{2}^{\left(1\right)} & \cdots & d_{n-3}^{\left(1\right)} & d_{n-2}^{\left(1\right)} & d_{n-1}^{\left(1\right)}\\
    d_{1} & d_{2} & d_{3} & \cdots & d_{n-2} & d_{n-1} & d_{n}\\
    d_{2}^{\left(-1\right)} & d_{3}^{\left(-1\right)} & d_{4}^{\left(-1\right)} & \cdots & d_{n-1}^{\left(-1\right)} & d_{n}^{\left(-1\right)} & 0\\
    d_{3}^{\left(-2\right)} & d_{4}^{\left(-2\right)} & d_{5}^{\left(-2\right)} & \cdots & d_{n}^{\left(-2\right)} & 0 & 0
    \end{matrix}\right)

The numbering here is a bit confusing, but in the column-wise storage, all
entries written in one column were in the same column in the original matrix.


Solving the system using pentapy
--------------------------------

To solve the system you can either provide :math:`M` as a full matrix or as
a flattened matrix in row-wise resp. col-wise flattened form.

If M is a full matrix, you call the following:

.. code-block:: python

    import pentapy as pp

    M = ...  # your matrix
    Y = ...  # your right hand side

    X = pp.solve(M, Y)

If M is flattend in row-wise order you have to set the keyword argument ``is_flat=True``:

.. code-block:: python

    import pentapy as pp

    M = ...  # your flattened matrix
    Y = ...  # your right hand side

    X = pp.solve(M, Y, is_flat=True)

If you got a col-wise flattend matrix you have to set ``index_row_wise=False``:

.. code-block:: python

    X = pp.solve(M, Y, is_flat=True, index_row_wise=False)


Tools
-----

pentapy provides some tools to convert a pentadiagonal matrix.

.. currentmodule:: pentapy.tools

.. autosummary::
   diag_indices
   shift_banded
   create_banded
   create_full


Gallery
=======



.. raw:: html

    <div class="sphx-glr-thumbnails">

.. thumbnail-parent-div-open

.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Here we create a random row wise flattened matrix M_flat and a random right hand side for the pentadiagonal equation system.">

.. only:: html

  .. image:: /examples/images/thumb/sphx_glr_01_solve_thumb.png
    :alt:

  :ref:`sphx_glr_examples_01_solve.py`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">1. Example: Solving</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Here we compare the outcome of the PTRANS-I and PTRANS-II algorithm for a random input.">

.. only:: html

  .. image:: /examples/images/thumb/sphx_glr_02_compare_thumb.png
    :alt:

  :ref:`sphx_glr_examples_02_compare.py`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">2. Example: Comparison</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Here we compare all algorithms for solving pentadiagonal systems provided by pentapy (except umf).">

.. only:: html

  .. image:: /examples/images/thumb/sphx_glr_03_perform_simple_thumb.png
    :alt:

  :ref:`sphx_glr_examples_03_perform_simple.py`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">3. Example: Performance for flattened matrices</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Here we compare all algorithms for solving pentadiagonal systems provided by pentapy (except umf) using a full quadratic matrix as input.">

.. only:: html

  .. image:: /examples/images/thumb/sphx_glr_04_perform_full_mat_thumb.png
    :alt:

  :ref:`sphx_glr_examples_04_perform_full_mat.py`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">4. Example: Performance for full matrices</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Here we demonstrate that the solver PTRANS-I can fail to solve a given system. A warning is given in that case and the output will be a nan-array.">

.. only:: html

  .. image:: /examples/images/thumb/sphx_glr_05_solve_error_thumb.png
    :alt:

  :ref:`sphx_glr_examples_05_solve_error.py`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">5. Example: Failing Corner Cases</div>
    </div>


.. thumbnail-parent-div-close

.. raw:: html

    </div>


.. toctree::
   :hidden:

   /examples/01_solve
   /examples/02_compare
   /examples/03_perform_simple
   /examples/04_perform_full_mat
   /examples/05_solve_error



.. only:: html

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

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