Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • kwant/python-mumps
  • kel85uk/mumpy
  • jcapriot/python-mumps
  • vlukes/python-mumps
4 results
Show changes
Commits on Source (99)
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1595 additions and 181 deletions
.gitattributes 0 → 100644
View file @ f727e905
# GitHub syntax highlighting
pixi.lock linguist-language=YAML linguist-generated=true
.gitignore
View file @ f727e905
# Project-specific
mumpy/_mumps.pxd
mumpy/mumps.pyx
*.DS_Store*
# Byte-compiled / optimized / DLL files
__pycache__/
......@@ -14,6 +12,7 @@ __pycache__/
.Python
env/
build/
build-install/
develop-eggs/
dist/
downloads/
......@@ -28,70 +27,20 @@ var/
.installed.cfg
*.egg
# PyInstaller
# Usually these files are written by a python script from a template
# before PyInstaller builds the exe, so as to inject date/other infos into it.
*.manifest
*.spec
# Installer logs
pip-log.txt
pip-delete-this-directory.txt
# Unit test / coverage reports
htmlcov/
.tox/
.coverage
.coverage.*
.cache
nosetests.xml
.pytest_cache
coverage.xml
*,cover
.hypothesis/
# Translations
*.mo
*.pot
# Django stuff:
*.log
local_settings.py
# Flask stuff:
instance/
.webassets-cache
# Scrapy stuff:
.scrapy
# Sphinx documentation
docs/_build/
# PyBuilder
target/
# IPython Notebook
.ipynb_checkpoints
# pyenv
.python-version
# celery beat schedule file
celerybeat-schedule
# dotenv
.env
# virtualenv
venv/
ENV/
# Spyder project settings
.spyderproject
# Rope project settings
.ropeproject
mumpy/*/*.c
mumpy/*.c
*.DS_Store*
.idea/
# pixi environments
.pixi
*.egg-info
.gitlab-ci.yml 0 → 100644
View file @ f727e905
image: ghcr.io/prefix-dev/pixi:0.29.0
variables:
FF_USE_FASTZIP: "true"
CACHE_COMPRESSION_LEVEL: "fastest"
pre-commit:
variables:
PRE_COMMIT_HOME: $CI_PROJECT_DIR/.pre-commit-cache
script:
- pixi run -e precommit git config --global --add safe.directory $CI_PROJECT_DIR
- pixi run -e precommit check
cache:
key: "$CI_JOB_NAME"
paths:
- .pre-commit-cache
- .pixi
test:
script:
- pixi run test
cache:
key: "$CI_JOB_NAME"
paths:
- .pixi
- build
build:
script:
- pixi run -e publish uv build --sdist
- pixi run -e publish twine check dist/*
artifacts:
paths:
- dist
expire_in: "7 days"
cache:
key: "$CI_JOB_NAME"
paths:
- .pixi
publish to test pypi:
needs:
- test
- build
- pre-commit
rules:
- if: '$CI_COMMIT_TAG =~ /^v[0-9]+\.[0-9]+\.[0-9]+.*\+test$/' # vX.Y.Z.post1+test
script:
- pixi run -e publish twine upload -u __token__ -p $PYPI_TEST_TOKEN --repository testpypi dist/*
cache:
key: "$CI_JOB_NAME"
paths:
- .pixi
publish to pypi:
needs:
- test
- build
- pre-commit
rules:
- if: '$CI_COMMIT_TAG =~ /^v[0-9]+\.[0-9]+\.[0-9]+[^+]*$/' # No +test
script:
- pixi run -e publish twine upload -u __token__ -p $PYPI_TOKEN dist/*
cache:
key: "$CI_JOB_NAME"
paths:
- .pixi
.pre-commit-config.yaml 0 → 100644
View file @ f727e905
repos:
- repo: https://github.com/charliermarsh/ruff-pre-commit
rev: "v0.6.4"
hooks:
- id: ruff
- id: ruff-format
- repo: https://github.com/PyCQA/doc8
rev: "v1.1.2"
hooks:
- id: doc8
- repo: https://github.com/codespell-project/codespell
rev: v2.3.0
hooks:
- id: codespell
stages: [commit]
additional_dependencies:
- tomli
# Prevent committing inline conflict markers
- repo: https://github.com/pre-commit/pre-commit-hooks
rev: v4.6.0
hooks:
- id: check-merge-conflict
args: [--assume-in-merge]
- id: check-yaml
- id: end-of-file-fixer
- id: trailing-whitespace
AUTHORS.md 0 → 100644
View file @ f727e905
# Python-MUMPS Authors
Below is a list of the contributors to Python-MUMPS
+ [Anton Akhmerov](<https://antonakhmerov.org>)
+ Christoph Groth
+ Kelvin Loh
+ [Bas Nijholt](<https://nijho.lt>)
+ [Joseph Weston](<https://joseph.weston.cloud>)
+ Michael Wimmer
+ Kostas Vilkelis
For a full list of contributors run
```bash
git log --pretty=format:"%an" | sort | uniq
```
LICENSE
View file @ f727e905
BSD 2-Clause License
Copyright (c) 2017, Anton Akhmerov, Christoph Groth, and Michael Wimmer,
Joseph Weston, and Bas Nijholt.
Copyright (c) 2018 Python-MUMPS Authors
All rights reserved.
......
README.md
View file @ f727e905
# mumpy
Python bindings for the MUMPS package MUMPS: a parallel sparse direct solver
# Python-MUMPS
# Installation
## On Linux
Either install `mumps` with:
```
sudo apt-get instal libmumps-scotch-dev
```
Python bindings for the [MUMPS](http://mumps-solver.org/): a parallel sparse direct solver.
## Both on Linux and OSX
or use the `conda` package on `conda-forge`
```
conda config --add channels conda-forge
conda install cython mumps numpy scipy libgfortran
pip install git+git@github.com:basnijholt/mumpy.git
```
and use the following `build.conf` (if on Linux, not needed for OSX)
## Scope
This package targets MUMPS packaged by conda-forge using Cython bindings. It
aims to provide a full wrapper of the MUMPS sequential API. Its primary target
OS is Linux.
Next steps include:
- Support for Windows and OSX
- Support for distributed (MPI) MUMPS
## Installation
`python-mumps` works with Python 3.10 and higher on Linux, Windows and Mac.
The recommended way to install `python-mumps` is using `mamba`/`conda`.
```bash
mamba install -c conda-forge python-mumps
```
[mumps]
include_dirs = $CONDA_PREFIX/include
library_dirs = $CONDA_PREFIX/lib
libraries = zmumps mumps_common pord metis esmumps scotch scotcherr mpiseq gfortran
extra_link_args = -Wl,-rpath=$CONDA_PREFIX/lib
```
# Usage
`python-mumps` can also be installed from PyPI, however this is a more involved procedure
that requires separately installing the MUMPS library and a C compiler.
## Usage example
The following example shows how Python-MUMPS can be used to implement sparse diagonalization
with Scipy.
```python
import scipy.sparse.linalg as sla
from scipy.sparse import identity
import mumpy.mumps as mumps
import mumps
def sparse_diag(matrix, k, sigma, **kwargs):
"""Call sla.eigsh with mumps support.
Please see scipy.sparse.linalg.eigsh for documentation.
Notes
-----
mumpy only works with complex numbers at the moment.
See scipy.sparse.linalg.eigsh for documentation.
"""
class LuInv(sla.LinearOperator):
def __init__(self, A):
inst = mumps.MUMPSContext()
inst = mumps.Context()
inst.analyze(A, ordering='pord')
inst.factor(A)
self.solve = inst.solve
......@@ -55,3 +56,31 @@ def sparse_diag(matrix, k, sigma, **kwargs):
opinv = LuInv(matrix - sigma * identity(matrix.shape[0]))
return sla.eigsh(matrix, k, sigma=sigma, OPinv=opinv, **kwargs)
```
## Development
### Pixi
`python-mumps` recommends [pixi](https://pixi.sh/).
After installing pixi, use
```bash
pixi run test -v # (Pytest arguments go after test)
```
This will also install the necessary dependencies.
### pre-commit
`python-mumps` uses [pre-commit](https://pre-commit.com/) to enforce code style. After installing it, run
```bash
pre-commit install
```
or if you want to use pre-commit provided by pixi, run
```bash
pixi run pre-commit install
```
meson.build 0 → 100644
View file @ f727e905
project(
'python-mumps',
'cython', 'c',
license: 'BSD-2-Clause',
license_files: 'LICENSE',
meson_version: '>= 1.1.0',
default_options: [
'buildtype=release',
'c_std=c99',
],
version: run_command('mumps/_version.py', check: true).stdout().strip(),
)
py = import('python').find_installation(pure: false)
cc = meson.get_compiler('c')
tempita = files('tempita.py')
fs = import('fs')
incdir_numpy = meson.get_external_property('numpy-include-dir', 'not-given')
if incdir_numpy == 'not-given'
incdir_numpy = run_command(py,
[
'-c',
'''import os
import numpy as np
try:
incdir = os.path.relpath(np.get_include())
except Exception:
incdir = np.get_include()
print(incdir)
'''
],
check: true
).stdout().strip()
else
_incdir_numpy_abs = incdir_numpy
endif
inc_np = include_directories(incdir_numpy)
numpy_nodepr_api = ['-DNPY_NO_DEPRECATED_API=NPY_1_9_API_VERSION']
np_dep = declare_dependency(include_directories: inc_np, compile_args: numpy_nodepr_api)
mumps_names = get_option('mumps_names')
mumps_libs = []
foreach mumps_name : mumps_names
if cc.find_library(mumps_name, required: false).found()
mumps_libs += cc.find_library(mumps_name)
else
# Fall back to the name without '_seq'.
mumps_libs += cc.find_library(mumps_name.substring(0, -4))
endif
endforeach
subdir('mumps')
meson.options 0 → 100644
View file @ f727e905
option(
'mumps_names',
type: 'array',
value: ['dmumps_seq', 'zmumps_seq', 'cmumps_seq', 'smumps_seq'],
description: 'MUMPS library names to link against. Default is: ones used in conda-forge.'
)
mumps/__init__.py 0 → 100644
View file @ f727e905
from .mumps import (
Context,
MUMPSError,
AnalysisStatistics,
FactorizationStatistics,
schur_complement,
nullspace,
possible_orderings,
orderings,
)
from ._version import __version__
__all__ = [
"Context",
"MUMPSError",
"AnalysisStatistics",
"FactorizationStatistics",
"schur_complement",
"nullspace",
"possible_orderings",
"orderings",
"__version__",
"test",
]
def test(verbose=True):
try:
from pytest import main
except ImportError:
raise RuntimeError("pytest is required to run the tests")
import os.path
return main(
[os.path.dirname(os.path.abspath(__file__)), "-s"] + (["-v"] if verbose else [])
)
test.__test__ = False
mumpy/_mumps.pxd.j2 → mumps/_mumps.pxd.in
View file @ f727e905
# Copyright 2011-2016 Anton Akhmerov, Christoph Groth, and Michael Wimmer
# Copyright 2017 Bas Nijholt.
# Copyright 2011-2016 Kwant authors.
# Copyright 2018 Python-MUMPS Authors.
#
# This file is part of mumpy. It is subject to the license terms in the file
# This file is part of Python-MUMPS. It is subject to the license terms in the file
# LICENSE found in the top-level directory of this distribution. A list of
# mumpy authors can be found in the file AUTHORS.md at the top-level
# directory of this distribution and at https://github.com/basnijholt/mumpy
# Python-MUMPS authors can be found in the file AUTHORS.md at the top-level
# directory of this distribution and at
# https://gitlab.kwant-project.org/kwant/python-mumps.
"""Cython declarations for the C interface to Mumps."""
cdef extern from "mumps_c_types.h":
cdef extern from "mumps_c_types.h" nogil:
ctypedef int MUMPS_INT
ctypedef float SMUMPS_REAL
......@@ -28,13 +29,13 @@ cdef extern from "mumps_c_types.h":
ctypedef double ZMUMPS_REAL
ctypedef mumps_double_complex ZMUMPS_COMPLEX
{% for p in 'SDCZ' %}
{{for p in 'SDCZ'}}
cdef extern from "{{p|lower}}mumps_c.h":
cdef extern from "{{p.lower()}}mumps_c.h" nogil:
ctypedef struct {{p}}MUMPS_STRUC_C:
MUMPS_INT sym, par, job
MUMPS_INT comm_fortran
MUMPS_INT icntl[40]
MUMPS_INT icntl[60]
{{p}}MUMPS_REAL cntl[15]
MUMPS_INT n
......@@ -61,6 +62,9 @@ cdef extern from "{{p|lower}}mumps_c.h":
MUMPS_INT *listvar_schur
{{p}}MUMPS_COMPLEX *schur
cdef void {{p|lower}}mumps_c({{p}}MUMPS_STRUC_C *)
MUMPS_INT lredrhs
{{p}}MUMPS_COMPLEX *redrhs
{% endfor %}
cdef void {{p.lower()}}mumps_c({{p}}MUMPS_STRUC_C *)
{{endfor}}
mumpy/mumps.pyx.j2 → mumps/_mumps.pyx.in
View file @ f727e905
# Copyright 2011-2016 Anton Akhmerov, Christoph Groth, and Michael Wimmer
# Copyright 2017 Joseph Weston and Bas Nijholt.
# Copyright 2011-2016 Kwant authors.
# Copyright 2018 Python-MUMPS Authors.
#
# This file is part of mumpy. It is subject to the license terms in the file
# This file is part of Python-MUMPS. It is subject to the license terms in the file
# LICENSE found in the top-level directory of this distribution. A list of
# mumpy authors can be found in the file AUTHORS.md at the top-level
# directory of this distribution and at https://github.com/basnijholt/mumpy.
# Python-MUMPS authors can be found in the file AUTHORS.md at the top-level
# directory of this distribution and at
# https://gitlab.kwant-project.org/kwant/python-mumps.
"""Low-level interface to Mumps.
Thin wrappers around the Mumps C API. At this stage there
......@@ -12,13 +13,21 @@ are still separate interfaces for real vs. complex and single
vs. double precision.
End users should not use this interface directly unless they
need finer control that that provided by the 'mumpy' interface.
need finer control that that provided by the 'mumps' interface.
"""
cimport numpy as np
import numpy as np
from . cimport _mumps as mumps
from .fortran_helpers import assert_fortran_matvec, assert_fortran_mat
from mumps cimport _mumps as mumps
from mumps.fortran_helpers import assert_fortran_matvec, assert_fortran_mat
from cpython.pythread cimport (
PyThread_type_lock,
PyThread_allocate_lock,
PyThread_acquire_lock,
PyThread_release_lock,
PyThread_free_lock
)
int_dtype = np.int32
......@@ -53,9 +62,9 @@ cdef make_mumps_int_array(mumps.MUMPS_INT *array):
wrapper.array = array
return wrapper
{% for p in 'SDCZ' %}
{{for p in 'SDCZ'}}
cdef class {{p|lower}}mumps_real_array:
cdef class {{p.lower()}}mumps_real_array:
cdef mumps.{{p}}MUMPS_REAL *array
def __init__(self):
......@@ -68,52 +77,65 @@ cdef class {{p|lower}}mumps_real_array:
self.array[key - 1] = value
cdef make_{{p|lower}}mumps_real_array(mumps.{{p}}MUMPS_REAL *array):
wrapper = {{p|lower}}mumps_real_array()
cdef make_{{p.lower()}}mumps_real_array(mumps.{{p}}MUMPS_REAL *array):
wrapper = {{p.lower()}}mumps_real_array()
wrapper.array = array
return wrapper
{% endfor %}
{{endfor}}
#############################################################
{% for p in 'SDCZ' %}
{{for p in 'SDCZ'}}
cdef class {{p|lower}}mumps:
cdef class {{p.lower()}}mumps:
cdef mumps.{{p}}MUMPS_STRUC_C params
cdef public mumps_int_array icntl
cdef public {{p|lower}}mumps_real_array cntl
cdef public {{p.lower()}}mumps_real_array cntl
cdef public mumps_int_array info
cdef public mumps_int_array infog
cdef public {{p|lower}}mumps_real_array rinfo
cdef public {{p|lower}}mumps_real_array rinfog
cdef public {{p.lower()}}mumps_real_array rinfo
cdef public {{p.lower()}}mumps_real_array rinfog
cdef PyThread_type_lock lock
def __cinit__(self, verbose=False, sym=0):
self.lock = PyThread_allocate_lock()
def __init__(self, verbose=False, sym=0):
self.params.job = -1
self.params.sym = sym
self.params.par = 1
self.params.comm_fortran = -987654
self.params.comm_fortran = -987654 # MPI_COMM_WORLD, according to docs
mumps.{{p|lower}}mumps_c(&self.params)
with nogil:
mumps.{{p.lower()}}mumps_c(&self.params)
self.icntl = make_mumps_int_array(self.params.icntl)
self.cntl = make_{{p|lower}}mumps_real_array(self.params.cntl)
self.cntl = make_{{p.lower()}}mumps_real_array(self.params.cntl)
self.info = make_mumps_int_array(self.params.info)
self.infog = make_mumps_int_array(self.params.infog)
self.rinfo = make_{{p|lower}}mumps_real_array(self.params.rinfo)
self.rinfog = make_{{p|lower}}mumps_real_array(self.params.rinfog)
self.rinfo = make_{{p.lower()}}mumps_real_array(self.params.rinfo)
self.rinfog = make_{{p.lower()}}mumps_real_array(self.params.rinfog)
def __init__(self, verbose=False, sym=0):
# no diagnostic output (MUMPS is very verbose normally)
if not verbose:
self.icntl[1] = 0
self.icntl[3] = 0
def __dealloc__(self):
self.params.job = -2
mumps.{{p|lower}}mumps_c(&self.params)
if self.lock: # guard to ensure __dealloc__ operations are not called multiple times
self.params.job = -2
with nogil:
mumps.{{p.lower()}}mumps_c(&self.params)
PyThread_free_lock(self.lock)
self.lock = NULL
def call(self):
mumps.{{p|lower}}mumps_c(&self.params)
with nogil:
PyThread_acquire_lock(self.lock, 1)
mumps.{{p.lower()}}mumps_c(&self.params)
PyThread_release_lock(self.lock)
def _set_job(self, value):
self.params.job = value
......@@ -131,7 +153,7 @@ cdef class {{p|lower}}mumps:
mumps.MUMPS_INT N,
np.ndarray[mumps.MUMPS_INT, ndim=1] i,
np.ndarray[mumps.MUMPS_INT, ndim=1] j,
np.ndarray[np.complex128_t, ndim=1] a):
np.ndarray[{{p}}MUMPS_COMPLEX_NPT, ndim=1] a):
self.params.n = N
self.params.nz = a.shape[0]
self.params.irn = <mumps.MUMPS_INT *>i.data
......@@ -141,8 +163,9 @@ cdef class {{p|lower}}mumps:
def set_dense_rhs(self, np.ndarray rhs):
assert_fortran_matvec(rhs)
if rhs.dtype != np.complex128:
raise ValueError("numpy array must be of dtype complex128!")
if rhs.dtype != {{p}}MUMPS_COMPLEX_NP:
raise ValueError("numpy array must be of dtype {}".format({{p}}MUMPS_COMPLEX_NP))
if rhs.ndim == 1:
self.params.nrhs = 1
......@@ -154,7 +177,7 @@ cdef class {{p|lower}}mumps:
def set_sparse_rhs(self,
np.ndarray[mumps.MUMPS_INT, ndim=1] col_ptr,
np.ndarray[mumps.MUMPS_INT, ndim=1] row_ind,
np.ndarray[np.complex128_t, ndim=1] data):
np.ndarray[{{p}}MUMPS_COMPLEX_NPT, ndim=1] data):
if row_ind.shape[0] != data.shape[0]:
raise ValueError("Number of entries in row index and value "
......@@ -167,7 +190,7 @@ cdef class {{p|lower}}mumps:
self.params.irhs_ptr = <mumps.MUMPS_INT *>col_ptr.data
def set_schur(self,
np.ndarray[np.complex128_t, ndim=2, mode='c'] schur,
np.ndarray[{{p}}MUMPS_COMPLEX_NPT, ndim=2, mode='c'] schur,
np.ndarray[mumps.MUMPS_INT, ndim=1] schur_vars):
if schur.shape[0] != schur.shape[1]:
......@@ -180,4 +203,10 @@ cdef class {{p|lower}}mumps:
self.params.schur = <mumps.{{p}}MUMPS_COMPLEX *>schur.data
self.params.listvar_schur = <mumps.MUMPS_INT *>schur_vars.data
{% endfor %}
def set_schur_rhs(self,
np.ndarray[{{p}}MUMPS_COMPLEX_NPT, ndim=1, mode='c'] rhs):
self.params.lredrhs = rhs.shape[0]
self.params.redrhs = <mumps.{{p}}MUMPS_COMPLEX *>rhs.data
{{endfor}}
mumps/_version.py 0 → 100755
View file @ f727e905
#!/usr/bin/env python
try:
# setuptools_scm is not a runtime dependency. This is used when pip
# installing from a git repository.
from setuptools_scm import get_version
__version__ = get_version(root="..", relative_to=__file__)
except (ImportError, LookupError):
from importlib.metadata import version, PackageNotFoundError
try:
__version__ = version("python-mumps")
except PackageNotFoundError:
__version__ = "0.0.0+unknown"
if __name__ == "__main__":
print(__version__)
mumpy/fortran_helpers.py → mumps/fortran_helpers.py
View file @ f727e905
# Copyright 2011-2016 Anton Akhmerov, Christoph Groth, Michael Wimmer
# Copyright 2017 Joseph Weston and Bas Nijholt.
# Copyright 2018 Python-MUMPS Authors
#
# This file is part of mumpy. It is subject to the license terms in the file
# This file is part of Python-MUMPS. It is subject to the license terms in the file
# LICENSE found in the top-level directory of this distribution. A list of
# mumpy authors can be found in the file AUTHORS.md at the top-level
# directory of this distribution and at https://github.com/basnijholt/mumpy
# Python-MUMPS authors can be found in the file AUTHORS.md at the top-level
# directory of this distribution and at
# https://gitlab.kwant-project.org/kwant/python-mumps.
import numpy as np
......@@ -30,18 +31,17 @@ def prepare_for_fortran(overwrite, *args):
determine the proper data type.
`prepare_for_lapack` returns a character indicating the proper
data type in LAPACK style ('s', 'd', 'c', 'z') and a list of
data type in LAPACK style 'SDCZ' and a list of
properly converted arrays.
"""
# Make sure we have NumPy arrays
mats = [None]*len(args)
mats = [None] * len(args)
for i in range(len(args)):
if args[i] is not None:
arr = np.asanyarray(args[i])
if not np.issubdtype(arr.dtype, np.number):
raise ValueError("Argument cannot be interpreted "
"as a numeric array")
raise ValueError("Argument cannot be interpreted as a numeric array")
mats[i] = (arr, arr is not args[i] or overwrite)
else:
......@@ -53,23 +53,23 @@ def prepare_for_fortran(overwrite, *args):
dtype = np.common_type(*[arr for arr, ovwrt in mats if arr is not None])
if dtype == np.float32:
lapacktype = 's'
lapacktype = "s"
elif dtype == np.float64:
lapacktype = 'd'
lapacktype = "d"
elif dtype == np.complex64:
lapacktype = 'c'
lapacktype = "c"
elif dtype == np.complex128:
lapacktype = 'z'
lapacktype = "z"
else:
raise AssertionError("Unexpected data type from common_type")
ret = [ lapacktype ]
ret = [lapacktype]
for npmat, ovwrt in mats:
# Now make sure that the array is contiguous, and copy if necessary.
if npmat is not None:
if npmat.ndim == 2:
if not npmat.flags["F_CONTIGUOUS"]:
npmat = np.asfortranarray(npmat, dtype = dtype)
npmat = np.asfortranarray(npmat, dtype=dtype)
elif npmat.dtype != dtype:
npmat = npmat.astype(dtype)
elif not ovwrt:
......@@ -78,7 +78,7 @@ def prepare_for_fortran(overwrite, *args):
npmat = np.asfortranarray(npmat.copy())
elif npmat.ndim == 1:
if not npmat.flags["C_CONTIGUOUS"]:
npmat = np.ascontiguousarray(npmat, dtype = dtype)
npmat = np.ascontiguousarray(npmat, dtype=dtype)
elif npmat.dtype != dtype:
npmat = npmat.astype(dtype)
elif not ovwrt:
......@@ -97,8 +97,11 @@ def assert_fortran_mat(*mats):
# This is a workaround for a bug in NumPy version < 2.0,
# where 1x1 matrices do not have the F_Contiguous flag set correctly.
for mat in mats:
if (mat is not None and (mat.shape[0] > 1 or mat.shape[1] > 1) and
not mat.flags["F_CONTIGUOUS"]):
if (
mat is not None
and (mat.shape[0] > 1 or mat.shape[1] > 1)
and not mat.flags["F_CONTIGUOUS"]
):
raise ValueError("Input matrix must be Fortran contiguous")
......@@ -109,11 +112,10 @@ def assert_fortran_matvec(*arrays):
# This is a workaround for a bug in NumPy version < 2.0,
# where 1x1 matrices do not have the F_Contiguous flag set correctly.
for arr in arrays:
if not arr.ndim in (1, 2):
raise ValueError("Input must be either a vector "
"or a matrix.")
if (not arr.flags["F_CONTIGUOUS"] or
(arr.ndim == 2 and arr.shape[0] == 1 and arr.shape[1] == 1) ):
raise ValueError("Input must be a Fortran ordered "
"NumPy array")
if arr.ndim not in (1, 2):
raise ValueError("Input must be either a vector or a matrix.")
if not arr.flags["F_CONTIGUOUS"] or (
arr.ndim == 2 and arr.shape[0] == 1 and arr.shape[1] == 1
):
raise ValueError("Input must be a Fortran ordered NumPy array")
mumps/meson.build 0 → 100644
View file @ f727e905
templates = files('_mumps.pxd.in', '_mumps.pyx.in')
pyx_files = []
foreach template : templates
base_name = fs.name(template)
rendered_name = base_name.replace('.in', '')
pyx_file = custom_target(
rendered_name,
input: template,
output: rendered_name,
command: [py, tempita, '@INPUT@', '@OUTPUT@']
)
pyx_files += pyx_file
endforeach
py_files = files('__init__.py', 'mumps.py', 'fortran_helpers.py', '_version.py')
py.install_sources(py_files, subdir: 'mumps')
foreach py_file : py_files
fs.copyfile(py_file)
endforeach
py.extension_module('_mumps',
pyx_files,
dependencies: [np_dep, mumps_libs],
install: true,
subdir: 'mumps'
)
subdir('tests')
mumpy/mumpy.py → mumps/mumps.py
View file @ f727e905
This diff is collapsed.
mumpy/tests/__init__.py → mumps/tests/__init__.py
View file @ f727e905
File moved
mumpy/tests/_test_utils.py → mumps/tests/_test_utils.py
View file @ f727e905
# Copyright 2011-2016 Anton Akhmerov, Christoph Groth, and Michael Wimmer
# Copyright 2017 Joseph Weston and Bas Nijholt.
# Copyright 2011-2016 Kwant authors.
# Copyright 2018 Python-MUMPS Authors.
#
# This file is part of mumpy. It is subject to the license terms in the file
# This file is part of Python-MUMPS. It is subject to the license terms in the file
# LICENSE found in the top-level directory of this distribution. A list of
# mumpy authors can be found in the file AUTHORS.md at the top-level
# directory of this distribution and at https://github.com/basnijholt/mumpy
# Python-MUMPS authors can be found in the file AUTHORS.md at the top-level
# directory of this distribution and at
# https://gitlab.kwant-project.org/kwant/python-mumps.
import numpy as np
......@@ -26,7 +27,7 @@ class _Random:
self._x = (a * self._x + c) % m
return (float(self._x)/m-0.5)*2
return (float(self._x) / m - 0.5) * 2
def _randi(self):
# A very bad random number generator returning number between 0 and 20.
......@@ -41,25 +42,25 @@ class _Random:
return self._x % 21
def randmat(self, n, m, dtype):
mat = np.empty((n, m), dtype = dtype)
mat = np.empty((n, m), dtype=dtype)
if issubclass(dtype, np.complexfloating):
for i in range(n):
for j in range(m):
mat[i,j] = self._randf() + 1j * self._randf()
mat[i, j] = self._randf() + 1j * self._randf()
elif issubclass(dtype, np.floating):
for i in range(n):
for j in range(m):
mat[i,j] = self._randf()
mat[i, j] = self._randf()
elif issubclass(dtype, np.integer):
for i in range(n):
for j in range(m):
mat[i,j] = self._randi()
mat[i, j] = self._randi()
return mat
def randvec(self, n, dtype):
vec = np.empty(n, dtype = dtype)
vec = np.empty(n, dtype=dtype)
if issubclass(dtype, np.complexfloating):
for i in range(n):
......@@ -72,14 +73,3 @@ class _Random:
vec[i] = self._randi()
return vec
# Improved version of assert_arrays_almost_equal that uses the dtype to set the
# precision (The default precision of assert_arrays_almost_equal is sometimes
# too small for single-precision comparisions.)
def assert_array_almost_equal(dtype, a, b):
if dtype == np.float32 or dtype == np.complex64:
prec = 5
else:
prec = 10
np.testing.assert_array_almost_equal(a, b, decimal=prec)
mumps/tests/meson.build 0 → 100644
View file @ f727e905
py.install_sources(
files('__init__.py', '_test_utils.py', 'test_mumps.py'),
subdir: 'mumps/tests'
)
mumps/tests/test_mumps.py 0 → 100644
View file @ f727e905
# Copyright 2011-2016 Kwant authors.
# Copyright 2018 Python-MUMPS Authors.
#
# This file is part of Python-MUMPS. It is subject to the license terms in the file
# LICENSE found in the top-level directory of this distribution. A list of
# Python-MUMPS authors can be found in the file AUTHORS.md at the top-level
# directory of this distribution and at
# https://gitlab.kwant-project.org/kwant/python-mumps.
import pytest
import numpy as np
import scipy.sparse as sp
import scipy.linalg as la
from mumps import Context, MUMPSError, schur_complement, nullspace
from ._test_utils import _Random
# Decimal places of precision per datatype. These limits have been determined
# heuristically by inspecting the upper error bound reported by MUMPS for
# the random matrices used here.
precisions = {
np.float32: 1, # yes, really only 1 decimal place of precision!
np.float64: 10,
np.complex64: 1,
np.complex128: 10,
}
dtypes = list(precisions.keys())
def assert_array_almost_equal(dtype, a, b):
np.testing.assert_almost_equal(a, b, decimal=precisions[dtype])
def assert_almost_equal(dtype, a, b):
np.testing.assert_almost_equal(a, b, decimal=precisions[dtype])
@pytest.mark.parametrize("dtype", dtypes, ids=str)
@pytest.mark.parametrize("mat_size", [2, 10, 100], ids=str)
def test_lu_with_dense(dtype, mat_size):
rand = _Random()
a = rand.randmat(mat_size, mat_size, dtype)
bmat = rand.randmat(mat_size, mat_size, dtype)
bvec = rand.randvec(mat_size, dtype)
ctx = Context()
ctx.factor(sp.coo_matrix(a))
xvec = ctx.solve(bvec)
xmat = ctx.solve(bmat)
assert_array_almost_equal(dtype, np.dot(a, xvec), bvec)
assert_array_almost_equal(dtype, np.dot(a, xmat), bmat)
# now "sparse" right hand side
xvec = ctx.solve(sp.csc_matrix(bvec.reshape(mat_size, 1)))
xmat = ctx.solve(sp.csc_matrix(bmat))
assert_array_almost_equal(dtype, np.dot(a, xvec), bvec.reshape(mat_size, 1))
assert_array_almost_equal(dtype, np.dot(a, xmat), bmat)
@pytest.mark.parametrize("dtype", dtypes, ids=str)
@pytest.mark.parametrize("mat_size", [5, 10], ids=str)
def test_schur_complement_with_dense(dtype, mat_size):
rand = _Random()
a = rand.randmat(mat_size, mat_size, dtype)
s = schur_complement(sp.coo_matrix(a), list(range(3)))
assert_array_almost_equal(dtype, np.linalg.inv(s), np.linalg.inv(a)[:3, :3])
@pytest.mark.parametrize("dtype", dtypes, ids=str)
@pytest.mark.parametrize("mat_size", [5, 10], ids=str)
@pytest.mark.parametrize("symmetric_matrix", [True, False], ids=str)
def test_schur_complement_solution(dtype, mat_size, symmetric_matrix):
rand = _Random()
a = rand.randmat(mat_size, mat_size, dtype)
if symmetric_matrix:
a = a + a.T
bvec = rand.randvec(mat_size, dtype)
ctx = Context()
ctx.set_matrix(a, symmetric=symmetric_matrix)
ctx.schur(range(3))
xvec = ctx.solve_schur(bvec)
assert_array_almost_equal(dtype, a @ xvec, bvec)
@pytest.mark.parametrize("dtype", dtypes, ids=str)
def test_factor_error(dtype):
"""Test that an error is raised if factor is asked to reuse missing analysis."""
a = sp.identity(10, dtype=dtype)
with pytest.raises(ValueError):
Context().factor(a, reuse_analysis=True)
@pytest.mark.parametrize("dtype", dtypes, ids=str)
def test_error_minus_19(dtype):
"""Test if MUMPSError -19 is properly caught by increasing memory"""
ctx = Context()
a = sp.eye(5000, dtype=dtype)
# Create the context so we can modify it for factorization
ctx.analyze(a)
# We ensure that this first call creates a -19 error by
# allocating only 1 MB for factorization
ctx.mumps_instance.icntl[23] = 1 # Memory upper bound set to 1MB
ctx.mumps_instance.icntl[14] = 1 # Initial memory relaxation to 1%
ctx.mumps_instance.job = 2
ctx.mumps_instance.call()
# ensure that we really don't have enough memory
assert ctx.mumps_instance.infog[1] == -19
# This call should not raise any errors as
# it would successfully allocate memory
Context().factor(a)
@pytest.mark.parametrize("dtype", dtypes, ids=str)
@pytest.mark.parametrize("mat_size", [20, 50], ids=str)
@pytest.mark.parametrize("symmetric_matrix", [True, False], ids=str)
def test_nullspace(dtype, mat_size, symmetric_matrix):
"""Tests the nullspace wrapper by creating a rank deficient matrix
with known deficiency
This test is to check whether the wrapper works. The algorithm
to obtain the nullspace basis within MUMPS is still rather
unclear with the use of these null pivot threshold. Suppose
that we have to live with this uncertainty, we at least create
a test case where we know for sure the right nullspace
dimensions and check if the subroutine has successfully
calculated the nullspace basis.
"""
if not symmetric_matrix:
test = sp.coo_matrix([[1, 2], [2, 4]])
try:
nullspace(test, False)
except MUMPSError as error:
if error.error == -37:
pytest.skip("Installed MUMPS does not support unsymmetric nullspace")
else:
raise
# We ensure the dimension of the nullspace row vectors by creating a random
# column vector, a of size (mat_size x 1)
rand = _Random()
a = sp.coo_matrix(rand.randmat(mat_size, 1, dtype))
# We then create the symmetric matrix by making this operation:
# A = a * a^T, where A has size (mat_size x mat_size). Also because the
# matrix is created by an outer product of two column vectors, the rank of
# A is always 1. Hence, we know the dimension of the nullspace row vectors
# is always mat_size - 1. Also, since the two column vectors are the same,
# A is then a symmetric matrix.
if symmetric_matrix:
b = a.copy()
else:
b = sp.coo_matrix(rand.randmat(mat_size, 1, dtype))
A = a * b.transpose()
A_original = a * b.transpose()
# Remove the lower triangular elements in order to comply with the
# symmetric solver of MUMPS
if symmetric_matrix:
for i in range(a.shape[0]):
for j in range(i):
A[i, j] = 0
A.eliminate_zeros()
# Start with a small pivot_threshold value
pivot_threshold = 1e-5
maximum_threshold = 1.0
# We repeatedly add the threshold until we get the full nullspace basis
while True:
X = nullspace(A, symmetric_matrix, pivot_threshold)
pivot_threshold += 2 * 1e-5
if (X.shape[1] == A.shape[0] - 1) or (pivot_threshold > maximum_threshold):
break
# We know that the row vector dimension of the returned nullspace
# has to be mat_size-1
assert X.shape[1] == A_original.shape[0] - 1
# Check that indeed the nullspace does indeed give 0 when multiplied
# with the input matrix
assert_almost_equal(dtype, la.norm(A_original * X), 0.0)
@pytest.mark.parametrize("dtype", dtypes, ids=str)
def test_one_by_one(dtype):
"""Test a 1x1 matrix.
This is a regression test for
https://gitlab.kwant-project.org/kwant/python-mumps/-/issues/18
"""
ctx = Context()
ctx.factor(sp.eye(1, dtype=dtype))
assert_almost_equal(dtype, ctx.solve(np.array([1], dtype=dtype)), 1)
@pytest.mark.parametrize("dtype", dtypes, ids=str)
def test_zero_size_rhs(dtype):
"""Test that a 0xn rhs can be solved."""
a = np.random.randn(10, 10).astype(dtype)
ctx = Context()
ctx.factor(a)
rhs = np.zeros((10, 0), dtype=dtype)
assert_almost_equal(dtype, ctx.solve(rhs), rhs)
@pytest.mark.parametrize("dtype", dtypes, ids=str)
def test_symmetric_matrix(dtype):
"""Test that a symmetric matrix can be solved."""
n = 10
a = np.random.randn(n, n).astype(dtype)
if np.iscomplexobj(a):
a += (1j * np.random.randn(n, n)).astype(dtype)
a += a.T
ctx = Context()
ctx.set_matrix(a, symmetric=True)
ctx.factor()
rhs = np.random.randn(n, 1).astype(dtype)
assert_almost_equal(dtype, ctx.solve(rhs), la.solve(a, rhs))