Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[EXTENDED] PR 26 Parallelized multiple right-hand side support, fully typed tools #27

Open
wants to merge 62 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
62 commits
Select commit Hold shift + click to select a range
96f67f6
pkg: made `.venv` in `.gitignore` more general
MothNik May 25, 2024
3deb05c
wip: [11] refactored algorithm 1 to handle multiple right-hand sides
MothNik Jun 7, 2024
328d9c9
feat: [11] added doctest runs to `pytest`
MothNik Jun 7, 2024
90bdeb1
tests: [11] replaced single tools test by scalable parametrized tools…
MothNik Jun 7, 2024
5080641
tests: [11] removed version from coverage
MothNik Jun 7, 2024
2f57609
lint: [11] fixed block comment lint error
MothNik Jun 7, 2024
8db6ba4
style: [11] improved headlines
MothNik Jun 8, 2024
1cdd525
wip: [11] formatted
MothNik Jun 8, 2024
ae656cf
test: [11] created conditioned banded matrix creator for testing purp…
MothNik Jun 8, 2024
26752f2
tests: [11] add doctested reference solver; made ill-conditioning mor…
MothNik Jun 8, 2024
a0c8849
feat: [11] made error messages informative
MothNik Jun 8, 2024
ee6943a
fix: [23] disabled cdivision to fix error handling on Python side
MothNik Jun 8, 2024
0284843
feat/refactor: [11] enabled multipe right-hand sides for solver I; im…
MothNik Jun 8, 2024
c32bec2
tests: [11] added shape check to doctest
MothNik Jun 8, 2024
e7d92c8
tests: [11] added extensive parametrized tests for solver I that also…
MothNik Jun 8, 2024
647012f
tests: [11] added more intermediate sizes for tests
MothNik Jun 8, 2024
8ab514f
feat/fix: [11] added cython annotations to build process; fixed wrong…
MothNik Jun 8, 2024
5086959
refactor/doc: fixed missing variable declarations; added clarifying c…
MothNik Jun 8, 2024
f062b88
style: [11] made - signs better readable
MothNik Jun 8, 2024
15c787e
doc/refactor: [11] improved docs and comments of algorithm I; removed…
MothNik Jun 8, 2024
ef7ec51
wip: [11] restructured chaotic alias comparisons
MothNik Jun 8, 2024
94d110c
feat: [11] finalized multiple right-hand side support (serial) on Cyt…
MothNik Jun 8, 2024
ba5945f
tests: [11] unified test for both pentapy algorithms, made tests a lo…
MothNik Jun 8, 2024
026e8ea
tests: [11] removed superficial unittest
MothNik Jun 8, 2024
f6fa4cb
misc: [11] fixed typos
MothNik Jun 8, 2024
9b3e122
style/refactor: [11] made core code readable to humans and not only a…
MothNik Jun 8, 2024
92abbbd
feat/refactor: [11] ensured all solvers behave the same in terms of o…
MothNik Jun 8, 2024
fb0cf2a
fix: [11] fixed coverage typo
MothNik Jun 8, 2024
c2e8032
lint: [11] linted cython files
MothNik Jun 8, 2024
b4190f8
package: [11] made requirements file-based and dynamic to make facili…
MothNik Jun 8, 2024
4c1e3a0
feat: [11] updated chagelog
MothNik Jun 8, 2024
b4f79e5
fix: [11] fixed wrong coverage exclude
MothNik Jun 8, 2024
a2edcb4
doc: [11] improved wording for preserving shape
MothNik Jun 8, 2024
007bc5a
fix: [11] fixed changelog
MothNik Jun 8, 2024
997d372
feat: [11] enable multi-threaded parallelism for multiple right-hand …
MothNik Jun 9, 2024
0cf9890
tests: [11] renamed tests for serial mode
MothNik Jun 9, 2024
9040bb7
test: [11] reduced load on tests; transitioned to template-based appr…
MothNik Jun 9, 2024
a0b3388
doc/style: [11] fully typed `tools`; improved documentation
MothNik Jun 9, 2024
b93102b
package: [11] included build information for parallelized solvers
MothNik Jun 9, 2024
3b015dc
tests: [11] finalised parallel solver tests
MothNik Jun 9, 2024
5e86a34
tests/fix: [11] fixed inter-os-incompatibility of doctests
MothNik Jun 9, 2024
712b0c7
test/fix: [11] ? really fixed the inter-os-problems ?
MothNik Jun 9, 2024
6fd57d4
feat: [11] updated changelog
MothNik Jun 9, 2024
88e30e6
fix/doc: [11] reverted change that caused overwrite of `mat`; added a…
MothNik Jun 9, 2024
b4ce4a2
docs: [11] removed doubled defaults from docstrings
MothNik Jun 9, 2024
e71ff9d
docs: [11] fixed docstring and spelling inconsistencies
MothNik Jun 9, 2024
aef7d0f
feat: [11] enabled future possible validation of the quality of the s…
MothNik Jun 10, 2024
86ebf88
refactor: [11] made internal `workers`-handling safer
MothNik Jun 10, 2024
b5f8b95
docs: [11] updated outdated `README`
MothNik Jun 10, 2024
87264d7
refactor: [11] removed `validate` and leveraged `info`; reduced test …
MothNik Jun 10, 2024
8db881f
doc: [11] updated outdated index
MothNik Jun 10, 2024
37c2fbc
refactor: [11] made Cython structure smarter for algorithm I
MothNik Jun 11, 2024
896b1cc
wip: [11] made use of Enums for checks; prepared for unifying solver …
MothNik Jun 11, 2024
2a6daf4
refactor/tests: [11] fully unified Cython implementations of the solv…
MothNik Jun 11, 2024
549422d
refactor/tests: [11] split up cluttered ptrans-solving into dedicated…
MothNik Jun 11, 2024
a31fe71
fix: [11] fixed broken coverage of unknown error (was not skipped by …
MothNik Jun 11, 2024
ac5759b
refactor: [11] made mu occupy the central column of the factorized ma…
MothNik Jun 12, 2024
0945501
fix:
MothNik Jul 21, 2024
d1433f9
refactor:
MothNik Jul 23, 2024
edb4f33
refactor:
MothNik Jul 23, 2024
5b42428
doc:
MothNik Jul 23, 2024
743549b
refactor:
MothNik Jul 23, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
42 changes: 23 additions & 19 deletions src/pentapy/core.py
Original file line number Diff line number Diff line change
Expand Up @@ -169,29 +169,33 @@ def solve(
if single_rhs:
rhs = rhs[:, np.newaxis]

try:
solver_func = (
psolver.penta_solver1
if solver_inter == pmodels.PentaSolverAliases.PTRANS_I
else psolver.penta_solver2
)
# the respective solver is chosen ...
solver_func = (
psolver.penta_solver1
if solver_inter == pmodels.PentaSolverAliases.PTRANS_I
else psolver.penta_solver2
)

# if there was only a 1D right-hand side, the result has to be flattened
sol, info = solver_func( # NOTE: info is for potential future validation
np.ascontiguousarray(mat_flat),
np.ascontiguousarray(rhs),
workers,
False, # NOTE: this can enable validation in the future
)
# ... and the solver is called
sol, info = solver_func(
MothNik marked this conversation as resolved.
Show resolved Hide resolved
np.ascontiguousarray(mat_flat),
np.ascontiguousarray(rhs),
workers,
)

if single_rhs:
sol = sol.ravel()
# in case of failure, the solver will return NaNs and issue a warning
if info > 0:
warnings.warn(
f"pentapy: {solver_inter.name} solver encountered singular matrix at "
f"row index {info - 1}. Returning NaNs."
)
sol = np.full(shape=rhs_og_shape, fill_value=np.nan)

return sol
# in case of success, the solution can be returned (reshaped if necessary)
if single_rhs:
sol = sol.ravel()

except ZeroDivisionError:
warnings.warn("pentapy: PTRANS-I not suitable for input-matrix.")
return np.full(shape=rhs_og_shape, fill_value=np.nan)
return sol

# Case 2: LAPACK's banded solver
elif solver_inter == pmodels.PentaSolverAliases.LAPACK:
Expand Down
2 changes: 0 additions & 2 deletions src/pentapy/solver.pxd
Original file line number Diff line number Diff line change
Expand Up @@ -3,14 +3,12 @@ cdef double[::, ::1] c_penta_solver1(
double[::, ::1] mat_flat,
double[::, ::1] rhs,
int workers,
bint validate,
int* info,
)

cdef double[::, ::1] c_penta_solver2(
double[::, ::1] mat_flat,
double[::, ::1] rhs,
int workers,
bint validate,
int* info,
)
64 changes: 49 additions & 15 deletions src/pentapy/solver.pyx
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
# cython: language_level=3, boundscheck=False, wraparound=False, cdivision=False
# cython: language_level=3, boundscheck=False, wraparound=False, cdivision=True

"""
This is a solver linear equation systems with a penta-diagonal matrix,
Expand All @@ -24,7 +24,6 @@ def penta_solver1(
double[::, ::1] mat_flat,
double[::, ::1] rhs,
int workers,
bint validate,
):

# NOTE: info is defined to be overwritten for possible future validations
Expand All @@ -36,7 +35,6 @@ def penta_solver1(
mat_flat,
rhs,
workers,
validate,
&info,
)
),
Expand All @@ -48,7 +46,6 @@ def penta_solver2(
double[::, ::1] mat_flat,
double[::, ::1] rhs,
int workers,
bint validate,
):

# NOTE: info is defined to be overwritten for possible future validations
Expand All @@ -60,7 +57,6 @@ def penta_solver2(
mat_flat,
rhs,
workers,
validate,
&info,
)
),
Expand All @@ -74,7 +70,6 @@ cdef double[::, ::1] c_penta_solver1(
double[::, ::1] mat_flat,
double[::, ::1] rhs,
int workers,
bint validate,
int* info,
):
"""
Expand All @@ -100,12 +95,16 @@ cdef double[::, ::1] c_penta_solver1(
# === Solving the system of equations ===

# first, the matrix is factorized
c_penta_factorize_algo_1(
info[0] = c_penta_factorize_algo_1(
&mat_flat[0, 0],
mat_n_cols,
&mat_factorized[0, 0],
)

# in case of a zero-division, the function exits early
if info[0] > 0:
return result

# then, all the right-hand sides are solved
for iter_col in prange(
rhs_n_cols,
Expand All @@ -120,11 +119,10 @@ cdef double[::, ::1] c_penta_solver1(
&result[0, iter_col],
)

info[0] = 0
return result


cdef void c_penta_factorize_algo_1(
cdef int c_penta_factorize_algo_1(
double* mat_flat,
int64_t mat_n_cols,
double* mat_factorized,
Expand Down Expand Up @@ -169,11 +167,19 @@ cdef void c_penta_factorize_algo_1(

# === Factorization ===

# NOTE: in the following mu is manually checked for zero-division to extract the
# proper value of ``info`` and exit early in case of failure;
# ``info`` is set to the row count where the error occured as for LAPACK ``pbtrf``

# First row
mu_i = mat_flat[mat_row_base_idx_2]
if mu_i == 0.0:
return 1

al_i_minus_1 = mat_flat[mat_row_base_idx_1] / mu_i
be_i_minus_1 = mat_flat[0] / mu_i


mat_factorized[0] = 0.0
mat_factorized[1] = mu_i
mat_factorized[2] = 0.0
Expand All @@ -183,6 +189,9 @@ cdef void c_penta_factorize_algo_1(
# Second row
ga_i = mat_flat[mat_row_base_idx_3 + 1]
mu_i = mat_flat[mat_row_base_idx_2 + 1] - al_i_minus_1 * ga_i
if mu_i == 0.0:
return 2

al_i = (mat_flat[mat_row_base_idx_1 + 1] - be_i_minus_1 * ga_i) / mu_i
be_i = mat_flat[1] / mu_i

Expand All @@ -198,6 +207,8 @@ cdef void c_penta_factorize_algo_1(
e_i = mat_flat[mat_row_base_idx_4 + iter_row]
ga_i = mat_flat[mat_row_base_idx_3 + iter_row] - al_i_minus_1 * e_i
mu_i = mat_flat[mat_row_base_idx_2 + iter_row] - be_i_minus_1 * e_i - al_i * ga_i
if mu_i == 0.0:
return iter_row + 1

al_i_plus_1 = (mat_flat[mat_row_base_idx_1 + iter_row] - be_i * ga_i) / mu_i
al_i_minus_1 = al_i
Expand All @@ -219,6 +230,9 @@ cdef void c_penta_factorize_algo_1(
e_i = mat_flat[mat_row_base_idx_4 + mat_n_cols - 2]
ga_i = mat_flat[mat_row_base_idx_3 + mat_n_cols - 2] - al_i_minus_1 * e_i
mu_i = mat_flat[mat_row_base_idx_2 + mat_n_cols - 2] - be_i_minus_1 * e_i - al_i * ga_i
if mu_i == 0.0:
return mat_n_cols - 1

al_i_plus_1 = (mat_flat[mat_row_base_idx_1 + mat_n_cols - 2] - be_i * ga_i) / mu_i

mat_factorized[fact_curr_base_idx] = e_i
Expand All @@ -231,14 +245,16 @@ cdef void c_penta_factorize_algo_1(
e_i = mat_flat[mat_row_base_idx_4 + mat_n_cols - 1]
ga_i = mat_flat[mat_row_base_idx_3 + mat_n_cols - 1] - al_i * e_i
mu_i = mat_flat[mat_row_base_idx_2 + mat_n_cols - 1] - be_i * e_i - al_i_plus_1 * ga_i
if mu_i == 0.0:
return mat_n_cols

mat_factorized[fact_curr_base_idx + 5] = e_i
mat_factorized[fact_curr_base_idx + 6] = mu_i
mat_factorized[fact_curr_base_idx + 7] = ga_i
mat_factorized[fact_curr_base_idx + 8] = 0.0
mat_factorized[fact_curr_base_idx + 9] = 0.0

return
return 0


cdef int c_solve_penta_from_factorize_algo_1(
Expand Down Expand Up @@ -335,7 +351,6 @@ cdef double[::, ::1] c_penta_solver2(
double[::, ::1] mat_flat,
double[::, ::1] rhs,
int workers,
bint validate,
int* info,
):
"""
Expand All @@ -361,11 +376,13 @@ cdef double[::, ::1] c_penta_solver2(
# === Solving the system of equations ===

# first, the matrix is factorized
c_penta_factorize_algo_2(
info[0] = c_penta_factorize_algo_2(
&mat_flat[0, 0],
mat_n_cols,
&mat_factorized[0, 0],
)
if info[0] > 0:
return result

# then, all the right-hand sides are solved
for iter_col in prange(
Expand All @@ -381,10 +398,9 @@ cdef double[::, ::1] c_penta_solver2(
&result[0, iter_col],
)

info[0] = 0
return result

cdef void c_penta_factorize_algo_2(
cdef int c_penta_factorize_algo_2(
double* mat_flat,
int64_t mat_n_cols,
double* mat_factorized,
Expand Down Expand Up @@ -430,9 +446,16 @@ cdef void c_penta_factorize_algo_2(

# === Factorization ===

# NOTE: in the following ps is manually checked for zero-division to extract the
# proper value of ``info`` and exit early in case of failure;
# ``info`` is set to the row count where the error occured as for LAPACK ``pbtrf``

# First row

ps_i = mat_flat[mat_row_base_idx_2 + mat_n_cols - 1]
if ps_i == 0.0:
return mat_n_cols

si_i_plus_1 = mat_flat[mat_row_base_idx_3 + mat_n_cols - 1] / ps_i
phi_i_plus_1 = mat_flat[mat_row_base_idx_4 + mat_n_cols - 1] / ps_i

Expand All @@ -446,6 +469,9 @@ cdef void c_penta_factorize_algo_2(
# Second row
rho_i = mat_flat[mat_row_base_idx_1 + mat_n_cols - 2]
ps_i = mat_flat[mat_row_base_idx_2 + mat_n_cols - 2] - si_i_plus_1 * rho_i
if ps_i == 0.0:
return mat_n_cols - 1

si_i = (mat_flat[mat_row_base_idx_3 + mat_n_cols - 2] - phi_i_plus_1 * rho_i) / ps_i
phi_i = mat_flat[mat_row_base_idx_4 + mat_n_cols - 2] / ps_i

Expand All @@ -461,6 +487,9 @@ cdef void c_penta_factorize_algo_2(
b_i = mat_flat[iter_row]
rho_i = mat_flat[mat_row_base_idx_1 + iter_row] - si_i_plus_1 * b_i
ps_i = mat_flat[mat_row_base_idx_2 + iter_row] - phi_i_plus_1 * b_i - si_i * rho_i
if ps_i == 0.0:
return iter_row + 1

si_i_minus_1 = (mat_flat[mat_row_base_idx_3 + iter_row] - phi_i * rho_i) / ps_i
si_i_plus_1 = si_i
si_i = si_i_minus_1
Expand All @@ -479,6 +508,9 @@ cdef void c_penta_factorize_algo_2(
b_i = mat_flat[1]
rho_i = mat_flat[mat_row_base_idx_1 + 1] - si_i_plus_1 * b_i
ps_i = mat_flat[mat_row_base_idx_2 + 1] - phi_i_plus_1 * b_i - si_i * rho_i
if ps_i == 0.0:
return 2

si_i_minus_1 = (mat_flat[mat_row_base_idx_3 + 1] - phi_i * rho_i) / ps_i
si_i_plus_1 = si_i
si_i = si_i_minus_1
Expand All @@ -493,14 +525,16 @@ cdef void c_penta_factorize_algo_2(
b_i = mat_flat[0]
rho_i = mat_flat[mat_row_base_idx_1 + 0] - si_i_plus_1 * b_i
ps_i = mat_flat[mat_row_base_idx_2 + 0] - phi_i * b_i - si_i * rho_i
if ps_i == 0.0:
return 1

mat_factorized[4] = b_i
mat_factorized[3] = rho_i
mat_factorized[2] = ps_i
mat_factorized[1] = 0.0
mat_factorized[0] = 0.0

return
return 0


cdef int c_solve_penta_from_factorize_algo_2(
Expand Down
6 changes: 1 addition & 5 deletions tests/templates.py
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@
# === Constants ===

SEED = 19_031_977
REF_WARNING_CONTENT = "not suitable for input-matrix."
REF_WARNING_CONTENT = "singular matrix at row index"
N_ROWS = [
3, # important edge case
4, # important edge case
Expand All @@ -27,12 +27,8 @@
51, # odd
100, # ...
101,
500,
501,
1_000,
1_001,
5_000,
5_001,
]
SOLVER_ALIASES_PTRANS_I = [1, "1", "pTrAnS-I"]
SOLVER_ALIASES_PTRANS_II = [2, "2", "pTrAnS-Ii"]
Expand Down