Contents

pentapy Quickstart

pentapy is a toolbox to deal with pentadiagonal matrices in Python and solve the corresponding linear equation systems.

Installation

The package can be installed via pip. On Windows you can install WinPython to get Python and pip running.

pip install pentapy

There are pre-built wheels for Linux, MacOS and Windows for most Python versions.

To get the scipy solvers running, you have to install scipy or you can use the extra argument:

pip install pentapy[all]

Instead of “all” you can also typ “scipy” or “umfpack”.

References

The solver is based on the algorithms PTRANS-I and PTRANS-II presented by Askar et al. 2015.

Examples

Solving a pentadiagonal linear equation system

This is an example of how to solve a LES with a pentadiagonal matrix.

import numpy as np
import pentapy as pp

size = 1000
# create a flattened pentadiagonal matrix
M_flat = (np.random.random((5, size)) - 0.5) * 1e-5
V = np.random.random(size) * 1e5
# solve the LES with M_flat as row-wise flattened matrix
X = pp.solve(M_flat, V, is_flat=True)

# create the corresponding matrix for checking
M = pp.create_full(M_flat, col_wise=False)
# calculate the error
print(np.max(np.abs(np.dot(M, X) - V)))

This should give something like:

4.257890395820141e-08

Performance

In the following, a couple of solvers for pentadiagonal systems are compared:

• Solver 1: Standard linear algebra solver of Numpy np.linalg.solve (link)

• Solver 2: scipy.sparse.linalg.spsolve (link)

• Solver 3: Scipy banded solver [scipy.linalg.solve_banded](scipy.github.io/devdocs/generated/scipy.linalg.solve_banded.html)

• Solver 4: pentapy.solve with solver=1

• Solver 5: pentapy.solve with solver=2

The performance plot was created with perfplot (link).

Requirements

• Numpy >= 1.14.5

Optional

• SciPy

• scikit-umfpack

License

MIT

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: \(M\cdot X = Y\), where \(M\) is a quadratic \(n\times n\) matrix given by:

\[\begin{split}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)\end{split}\]

The aim is now, to solve this equation for \(X\).

Memory efficient storage

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

1. row-wise storage:

\[\begin{split}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)\end{split}\]

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:

\[\begin{split}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)\end{split}\]

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 \(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:

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:

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:

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

Tools

pentapy provides some tools to convert a pentadiagonal matrix.

diag_indices(n[, offset])	Get indices for the main or minor diagonals of a matrix.
shift_banded(mat[, up, low, col_to_row, copy])	Shift rows of a banded matrix.
create_banded(mat[, up, low, col_wise, dtype])	Create a banded matrix from a given quadratic Matrix.
create_full(mat[, up, low, col_wise])	Create a (n x n) Matrix from a given banded matrix.

Gallery

1. Example: Solving

1. Example: Solving

2. Example: Comparison

2. Example: Comparison

3. Example: Performance for flattened matrices

3. Example: Performance for flattened matrices

4. Example: Performance for full matrices

4. Example: Performance for full matrices

5. Example: Failing Corner Cases

5. Example: Failing Corner Cases

1. Example: Solving

Here we create a random row wise flattened matrix M_flat and a random right hand side for the pentadiagonal equation system.

After solving we calculate the absolute difference between the right hand side and the product of the matrix (which is transformed to a full quadratic one) and the solution of the system.

2.3836037144064903e-08

import numpy as np

import pentapy as pp

size = 1000
# create a flattened pentadiagonal matrix
M_flat = (np.random.random((5, size)) - 0.5) * 1e-5
V = np.random.random(size) * 1e5
# solve the LES with M_flat as row-wise flattened matrix
X = pp.solve(M_flat, V, is_flat=True)

# create the corresponding matrix for checking
M = pp.create_full(M_flat, col_wise=False)
# calculate the error
print(np.max(np.abs(np.dot(M, X) - V)))

2. Example: Comparison

Here we compare the outcome of the PTRANS-I and PTRANS-II algorithm for a random input.

rel. diff X1 X2:  1.1304241370505336e-14
max X1:  5226536599276.376
max X2:  5226536599276.326
max |M*X1 - V|:  1.837179297581315e-07
max |M*X2 - V|:  2.35304469242692e-08

import numpy as np

import pentapy as pp

size = 1000
# create a flattened pentadiagonal matrix
M_flat = (np.random.random((5, size)) - 0.5) * 1e-5
V = np.random.random(size) * 1e5
# compare the two solvers
X1 = pp.solve(M_flat, V, is_flat=True, solver=1)
X2 = pp.solve(M_flat, V, is_flat=True, solver=2)

# calculate the error
print("rel. diff X1 X2: ", np.max(np.abs(X1 - X2)) / np.max(np.abs(X1 + X2)))
print("max X1: ", np.max(np.abs(X1)))
print("max X2: ", np.max(np.abs(X2)))

M = pp.create_full(M_flat, col_wise=False)
# calculate the error
print("max |M*X1 - V|: ", np.max(np.abs(np.dot(M, X1) - V)))
print("max |M*X2 - V|: ", np.max(np.abs(np.dot(M, X2) - V)))

3. Example: Performance for flattened matrices

Here we compare all algorithms for solving pentadiagonal systems provided by pentapy (except umf).

To use this script you need to have the following packages installed:

• scipy

• perfplot

• matplotlib

Overall ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 100% 0:00:00
Kernels ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 100% 0:00:00

import numpy as np
import perfplot

from pentapy import solve, tools


def get_les(size):
    mat = (np.random.random((5, size)) - 0.5) * 1e-5
    V = np.array(np.random.random(size) * 1e5)
    return mat, V


def solve_1(in_val):
    """PTRANS-I"""
    mat, V = in_val
    return solve(mat, V, is_flat=True, index_row_wise=False, solver=1)


def solve_2(in_val):
    """PTRANS-II"""
    mat, V = in_val
    return solve(mat, V, is_flat=True, index_row_wise=False, solver=2)


def solve_3(in_val):
    mat, V = in_val
    return solve(mat, V, is_flat=True, index_row_wise=False, solver=3)


def solve_4(in_val):
    mat, V = in_val
    return solve(mat, V, is_flat=True, index_row_wise=False, solver=4)


def solve_5(in_val):
    mat, V = in_val
    M = tools.create_full(mat)
    return np.linalg.solve(M, V)


perfplot.show(
    setup=get_les,
    kernels=[solve_1, solve_2, solve_3, solve_4, solve_5],
    labels=[
        "PTRANS-I",
        "PTRANS-II",
        "scipy.linalg.solve_banded",
        "scipy.sparse.linalg.spsolve",
        "numpy.linalg.solve",
    ],
    n_range=[2**k for k in range(2, 13)],
    xlabel="Size [n]",
    logy=True,
)

4. Example: Performance for full matrices

Here we compare all algorithms for solving pentadiagonal systems provided by pentapy (except umf) using a full quadratic matrix as input.

To use this script you need to have the following packages installed:

• scipy

• perfplot

• matplotlib

Overall ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 100% 0:00:00
Kernels ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 100% 0:00:00

import numpy as np
import perfplot

from pentapy import solve, tools


def get_les(size):
    mat = (np.random.random((5, size)) - 0.5) * 1e-5
    V = np.array(np.random.random(size) * 1e5)
    M = tools.create_full(mat)
    return M, V


def solve_1(in_val):
    """PTRANS-I"""
    mat, V = in_val
    return solve(mat, V, is_flat=False, solver=1)


def solve_2(in_val):
    """PTRANS-II"""
    mat, V = in_val
    return solve(mat, V, is_flat=False, solver=2)


def solve_3(in_val):
    mat, V = in_val
    return solve(mat, V, is_flat=False, solver=3)


def solve_4(in_val):
    mat, V = in_val
    return solve(mat, V, is_flat=False, solver=4)


def solve_5(in_val):
    mat, V = in_val
    return np.linalg.solve(mat, V)


perfplot.show(
    setup=get_les,
    kernels=[solve_1, solve_2, solve_3, solve_4, solve_5],
    labels=[
        "PTRANS-I",
        "PTRANS-II",
        "scipy.linalg.solve_banded",
        "scipy.sparse.linalg.spsolve",
        "numpy.linalg.solve",
    ],
    n_range=[2**k for k in range(2, 13)],
    xlabel="Size [n]",
    logy=True,
)

5. Example: Failing Corner Cases

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.

[nan nan nan nan]

import numpy as np

import pentapy as pp

# create a full pentadiagonal matrix
mat = np.array([[3, 2, 1, 0], [-3, -2, 7, 1], [3, 2, -1, 5], [0, 1, 2, 3]])
V = np.array([6, 3, 9, 6])
# solve the LES with mat as a qudratic input matrix
X = pp.solve(mat, V, is_flat=False, solver=1)
print(X)

pentapy API

Purpose

pentapy is a toolbox to deal with pentadiagonal matrizes in Python.

Subpackages

tools

The tools module of pentapy.

Solver

Solver for a pentadiagonal equations system.

solve(mat, rhs[, is_flat, index_row_wise, ...])

Solver for a pentadiagonal system.

Tools

The following tools are provided:

diag_indices(n[, offset])	Get indices for the main or minor diagonals of a matrix.
shift_banded(mat[, up, low, col_to_row, copy])	Shift rows of a banded matrix.
create_banded(mat[, up, low, col_wise, dtype])	Create a banded matrix from a given quadratic Matrix.
create_full(mat[, up, low, col_wise])	Create a (n x n) Matrix from a given banded matrix.

Changelog

All notable changes to pentapy will be documented in this file.

1.2.0 - 2023-04

See #19

Enhancements

• added support for python 3.10 and 3.11

• add wheels for arm64 systems

• created solver.pxd file to be able to cimport the solver module

• added a CITATION.bib file

Changes

• move to src/ based package structure

• dropped python 3.7 support

• move meta-data to pyproject.toml

• simplified documentation

Bugfixes

• determine correct version when installing from archive

1.1.2 - 2021-07

Changes

• new package structure with pyproject.toml (#15)

• Sphinx-Gallery for Examples

• Repository restructuring: use a single main branch

• use np.asarray in solve to speed up computation (#17)

1.1.1 - 2021-02

Enhancements

• Python 3.9 support

Changes

• GitHub Actions for CI

1.1.0 - 2020-03-22

Enhancements

• Python 3.8 support

Changes

• python only builds are no longer available

• Python 2.7 and 3.4 support dropped

1.0.3 - 2019-11-10

Enhancements

• the algorithms PTRANS-I and PTRANS-II now raise a warning when they can not solve the given system

• there are now switches to install scipy and umf solvers as extra requirements

Bugfixes

• multiple minor bugfixes

1.0.0 - 2019-09-18

Enhancements

• the second algorithm PTRANS-II from Askar et al. 2015 is now implemented and can be used by solver=2

• the package is now tested and a coverage is calculated

• there are now pre-built binaries for Python 3.7

• the documentation is now available under https://geostat-framework.readthedocs.io/projects/pentapy

Changes

• pentapy is now licensed under the MIT license

0.1.1 - 2019-03-08

Bugfixes

• MANIFEST.in was missing in the 0.1.0 version

0.1.0 - 2019-03-07

This is the first release of pentapy, a python toolbox for solving pentadiagonal linear equation systems. The solver is implemented in cython, which makes it really fast.