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 (2.7, 3.4-3.7).
If your system is not supported and you want to have the Cython routines of pentapy installed, you have to provide a c-compiler and run:
pip install numpy cython
pip install pentapy
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).
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: , where
is a quadratic
matrix given by:
The aim is now, to solve this equation for .
Memory efficient storage¶
To store a pentadiagonal matrix memory efficient, there are two options:
row-wise storage:
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.
column-wise storage:
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 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.
|
Get indices for the main or minor diagonals of a matrix. |
|
Shift rows of a banded matrix. |
|
Create a banded matrix from a given quadratic Matrix. |
|
Create a (n x n) Matrix from a given banded matrix. |
Gallery¶
Note
Click here to download the full example code
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.
Out:
6.700440280837938e-06
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)))
Total running time of the script: ( 0 minutes 0.006 seconds)
Note
Click here to download the full example code
2. Example: Comparison¶
Here we compare the outcome of the PTRANS-I and PTRANS-II algorithm for a random input.
Out:
rel. diff X1 X2: 5.27456543735397e-15
max X1: 8146405710639.143
max X2: 8146405710639.057
max |M*X1 - V|: 4.2506144382059574e-08
max |M*X2 - V|: 1.2479722499847412e-07
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)))
Total running time of the script: ( 0 minutes 0.006 seconds)
Note
Click here to download the full example code
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

Out:
Overall ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 100% 0:00:00
Kernels 0% 0:00:01
from __future__ import division, absolute_import, print_function
import numpy as np
from pentapy import solve
from pentapy import tools
import perfplot
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,
)
Total running time of the script: ( 0 minutes 31.200 seconds)
Note
Click here to download the full example code
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

Out:
Overall ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 100% 0:00:00
Kernels 0% 0:00:01
from __future__ import division, absolute_import, print_function
import numpy as np
from pentapy import solve
from pentapy import tools
import perfplot
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,
)
Total running time of the script: ( 0 minutes 30.177 seconds)
Note
Click here to download the full example code
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.
Out:
[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)
Total running time of the script: ( 0 minutes 0.001 seconds)
pentapy API¶
Purpose¶
pentapy is a toolbox to deal with pentadiagonal matrizes in Python.
Solver¶
Solver for a pentadiagonal equations system.
|
Solver for a pentadiagonal system. |
Tools¶
The following tools are provided:
|
Get indices for the main or minor diagonals of a matrix. |
|
Shift rows of a banded matrix. |
|
Create a banded matrix from a given quadratic Matrix. |
|
Create a (n x n) Matrix from a given banded matrix. |
pentapy.core¶
The core module of pentapy.
The following functions are provided
|
Solver for a pentadiagonal system. |
-
pentapy.core.
solve
(mat, rhs, is_flat=False, index_row_wise=True, solver=1)[source]¶ Solver for a pentadiagonal system.
The matrix can be given as a full n x n matrix or as a flattend one. The flattend matrix can be given in a row-wise flattend form:
[[Dup2[0] Dup2[1] Dup2[2] ... Dup2[N-2] 0 0 ] [Dup1[0] Dup1[1] Dup1[2] ... Dup1[N-2] Dup1[N-1] 0 ] [Diag[0] Diag[1] Diag[2] ... Diag[N-2] Diag[N-1] Diag[N] ] [0 Dlow1[1] Dlow1[2] ... Dlow1[N-2] Dlow1[N-1] Dlow1[N]] [0 0 Dlow2[2] ... Dlow2[N-2] Dlow2[N-2] Dlow2[N]]]
Or a column-wise flattend form:
[[0 0 Dup2[2] ... Dup2[N-2] Dup2[N-1] Dup2[N] ] [0 Dup1[1] Dup1[2] ... Dup1[N-2] Dup1[N-1] Dup1[N] ] [Diag[0] Diag[1] Diag[2] ... Diag[N-2] Diag[N-1] Diag[N] ] [Dlow1[0] Dlow1[1] Dlow1[2] ... Dlow1[N-2] Dlow1[N-1] 0 ] [Dlow2[0] Dlow2[1] Dlow2[2] ... Dlow2[N-2] 0 0 ]]
Dup1 and Dup2 are the first and second upper minor-diagonals and Dlow1 resp. Dlow2 are the lower ones. If you provide a column-wise flattend matrix, you have to set:
index_row_wise=False
- Parameters
mat (
numpy.ndarray
) – The Matrix or the flattened Version of the pentadiagonal matrix.rhs (
numpy.ndarray
) – The right hand side of the equation system.is_flat (
bool
, optional) – State if the matrix is already flattend. Default:False
index_row_wise (
bool
, optional) – State if the flattend matrix is row-wise flattend. Default:True
solver (
int
orstr
, optional) –Which solver should be used. The following are provided:
[1, "1", "PTRANS-I"]
: The PTRANS-I algorithm[2, "2", "PTRANS-II"]
: The PTRANS-II algorithm[3, "3", "lapack", "solve_banded"]
: scipy.linalg.solve_banded[4, "4", "spsolve"]
: The scipy sparse solver without umf_pack[5, "5", "spsolve_umf", "umf", "umf_pack"]
: The scipy sparse solver with umf_pack
Default:
1
- Returns
result – Solution of the equation system
- Return type
pentapy.tools¶
The tools module of pentapy.
The following functions are provided
|
Get indices for the main or minor diagonals of a matrix. |
|
Shift rows of a banded matrix. |
|
Create a banded matrix from a given quadratic Matrix. |
|
Create a (n x n) Matrix from a given banded matrix. |
-
pentapy.tools.
create_banded
(mat, up=2, low=2, col_wise=True, dtype=None)[source]¶ Create a banded matrix from a given quadratic Matrix.
The Matrix will to be returned as a flattend matrix. Either in a column-wise flattend form:
[[0 0 Dup2[2] ... Dup2[N-2] Dup2[N-1] Dup2[N] ] [0 Dup1[1] Dup1[2] ... Dup1[N-2] Dup1[N-1] Dup1[N] ] [Diag[0] Diag[1] Diag[2] ... Diag[N-2] Diag[N-1] Diag[N] ] [Dlow1[0] Dlow1[1] Dlow1[2] ... Dlow1[N-2] Dlow1[N-1] 0 ] [Dlow2[0] Dlow2[1] Dlow2[2] ... Dlow2[N-2] 0 0 ]]
Then use:
col_wise=True
Or in a row-wise flattend form:
[[Dup2[0] Dup2[1] Dup2[2] ... Dup2[N-2] 0 0 ] [Dup1[0] Dup1[1] Dup1[2] ... Dup1[N-2] Dup1[N-1] 0 ] [Diag[0] Diag[1] Diag[2] ... Diag[N-2] Diag[N-1] Diag[N] ] [0 Dlow1[1] Dlow1[2] ... Dlow1[N-2] Dlow1[N-1] Dlow1[N]] [0 0 Dlow2[2] ... Dlow2[N-2] Dlow2[N-2] Dlow2[N]]]
Then use:
col_wise=False
Dup1 and Dup2 or the first and second upper minor-diagonals and Dlow1 resp. Dlow2 are the lower ones. The number of upper and lower minor-diagonals can be altered.
- Parameters
mat (
numpy.ndarray
) – The full (n x n) Matrix.up (
int
) – The number of upper minor-diagonals. Default: 2low (
int
) – The number of lower minor-diagonals. Default: 2col_wise (
bool
, optional) – Use column-wise storage. If False, use row-wise storage. Default:True
- Returns
Bandend matrix
- Return type
-
pentapy.tools.
create_full
(mat, up=2, low=2, col_wise=True)[source]¶ Create a (n x n) Matrix from a given banded matrix.
The given Matrix has to be a flattend matrix. Either in a column-wise flattend form:
[[0 0 Dup2[2] ... Dup2[N-2] Dup2[N-1] Dup2[N] ] [0 Dup1[1] Dup1[2] ... Dup1[N-2] Dup1[N-1] Dup1[N] ] [Diag[0] Diag[1] Diag[2] ... Diag[N-2] Diag[N-1] Diag[N] ] [Dlow1[0] Dlow1[1] Dlow1[2] ... Dlow1[N-2] Dlow1[N-1] 0 ] [Dlow2[0] Dlow2[1] Dlow2[2] ... Dlow2[N-2] 0 0 ]]
Then use:
col_wise=True
Or in a row-wise flattend form:
[[Dup2[0] Dup2[1] Dup2[2] ... Dup2[N-2] 0 0 ] [Dup1[0] Dup1[1] Dup1[2] ... Dup1[N-2] Dup1[N-1] 0 ] [Diag[0] Diag[1] Diag[2] ... Diag[N-2] Diag[N-1] Diag[N] ] [0 Dlow1[1] Dlow1[2] ... Dlow1[N-2] Dlow1[N-1] Dlow1[N]] [0 0 Dlow2[2] ... Dlow2[N-2] Dlow2[N-2] Dlow2[N]]]
Then use:
col_wise=False
Dup1 and Dup2 or the first and second upper minor-diagonals and Dlow1 resp. Dlow2 are the lower ones. The number of upper and lower minor-diagonals can be altered.
- Parameters
mat (
numpy.ndarray
) – The flattened Matrix.up (
int
) – The number of upper minor-diagonals. Default: 2low (
int
) – The number of lower minor-diagonals. Default: 2col_wise (
bool
, optional) – Input is in column-wise storage. If False, use as row-wise storage. Default:True
- Returns
Full matrix.
- Return type
-
pentapy.tools.
diag_indices
(n, offset=0)[source]¶ Get indices for the main or minor diagonals of a matrix.
This returns a tuple of indices that can be used to access the main diagonal of an array a with
a.ndim == 2
dimensions and shape (n, n).- Parameters
- Returns
idx (
numpy.ndarray
) – row indicesidy (
numpy.ndarray
) – col indices
-
pentapy.tools.
shift_banded
(mat, up=2, low=2, col_to_row=True, copy=True)[source]¶ Shift rows of a banded matrix.
Either from column-wise to row-wise storage or vice versa.
The Matrix has to be given as a flattend matrix. Either in a column-wise flattend form:
[[0 0 Dup2[2] ... Dup2[N-2] Dup2[N-1] Dup2[N] ] [0 Dup1[1] Dup1[2] ... Dup1[N-2] Dup1[N-1] Dup1[N] ] [Diag[0] Diag[1] Diag[2] ... Diag[N-2] Diag[N-1] Diag[N] ] [Dlow1[0] Dlow1[1] Dlow1[2] ... Dlow1[N-2] Dlow1[N-1] 0 ] [Dlow2[0] Dlow2[1] Dlow2[2] ... Dlow2[N-2] 0 0 ]]
Then use:
col_to_row=True
Or in a row-wise flattend form:
[[Dup2[0] Dup2[1] Dup2[2] ... Dup2[N-2] 0 0 ] [Dup1[0] Dup1[1] Dup1[2] ... Dup1[N-2] Dup1[N-1] 0 ] [Diag[0] Diag[1] Diag[2] ... Diag[N-2] Diag[N-1] Diag[N] ] [0 Dlow1[1] Dlow1[2] ... Dlow1[N-2] Dlow1[N-1] Dlow1[N]] [0 0 Dlow2[2] ... Dlow2[N-2] Dlow2[N-2] Dlow2[N]]]
Then use:
col_to_row=False
Dup1 and Dup2 are the first and second upper minor-diagonals and Dlow1 resp. Dlow2 are the lower ones. The number of upper and lower minor-diagonals can be altered.
- Parameters
mat (
numpy.ndarray
) – The Matrix or the flattened Version of the pentadiagonal matrix.up (
int
) – The number of upper minor-diagonals. Default: 2low (
int
) – The number of lower minor-diagonals. Default: 2col_to_row (
bool
, optional) – Shift from column-wise to row-wise storage or vice versa. Default:True
copy (
bool
, optional) – Copy the input matrix or overwrite it. Default:True
- Returns
Shifted bandend matrix
- Return type
Changelog¶
All notable changes to pentapy will be documented in this file.
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
andPTRANS-II
now raise a warning when they can not solve the given systemthere 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 bysolver=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