r"""Module with an interface to the Pardiso sparse solver.
>>> p = Pardiso()
>>> indptr = np.array([0, 4, 7, 9, 11, 14, 16, 17, 18])
>>> indices = np.array([0, 2, 5, 6,
... 1, 2, 4,
... 2, 7,
... 3, 6,
... 4, 5, 6,
... 5, 7,
... 6,
... 7])
>>> data = np.array([7.0, 1.0, 2.0, 7.0,
... -4.0, 8.0, 2.0,
... 1.0, 5.0,
... 7.0, 9.0,
... 5.0, 1.0, 5.0,
... 0.0, 5.0,
... 11.0,
... 5.0])
>>> A = sp.sparse.csr_matrix((data, indices, indptr))
>>> b = np.ones(A.shape[0], dtype=A.dtype)
>>> p.solve(A=A, b=b, symmetric=True)
... # doctest: +ELLIPSIS, +NORMALIZE_WHITESPACE
array([-0.0427..., -0.0158... , 0.1105..., -0.1204..., 0.0260...,
-0.1224..., 0.2047..., 0.2118...])
An asymmetric example:
>>> indptr = np.array([0, 4, 7, 9, 11, 12, 15, 17, 20])
>>> indices = np.array([0, 2, 5, 6,
... 1, 2, 4,
... 2, 7,
... 3, 6,
... 1,
... 2, 5, 7,
... 1, 6,
... 2, 6, 7])
>>> data = np.array([7.0, 1.0, 2.0, 7.0,
... -4.0, 8.0, 2.0,
... 1.0, 5.0,
... 7.0, 9.0,
... -4.0,
... 7.0, 3.0, 8.0,
... 1.0, 11.0,
... -3.0, 2.0, 5.0])
>>> p.solve(A=A, b=b, symmetric=False)
... # doctest: +ELLIPSIS, +NORMALIZE_WHITESPACE
array([ 0. , -0.213..., 0. , 0.0259..., 0.07272...,
0.1818..., 0.0909..., 0.2 ])
"""
from __future__ import division
__all__ = ['PardisoError', 'Pardiso']
import numpy as np
import scipy as sp
import scipy.sparse
import scipy.weave
from mmf.objects import StateVars, process_vars
[docs]class PardisoError(Exception):
"""Error with Pardiso solver."""
[docs]class Pardiso(StateVars):
r"""Encapsulate the Pardiso solver.
Notes
-----
The solver expects the input for matrices in the
:class:`scipy.sparse.csr_matrix` format. If the matrices are
symmetric, only the upper triangular portion should be provided.
"""
_state_vars = [
('solver', 0,
r"""This scalar value defines the solver method that the
user would like to use:
`0` : Sparse direct solver
`1` : Multi-recursive iterative solver"""),
('threads', None,
r"""Number of threads to use. If `None`, then look at the
environment variable `OMP_NUM_THREADS` or default to
1."""),
('verbose', False,
r"""If `True`, then pardiso will print out a buch of
statistics."""),
]
process_vars()
[docs] def __init__(self, *v, **kw):
if self.threads is None:
self.threads = 1
[docs] def mtype(self, A,
structurally_symmetric=False,
symmetric=False,
hermitian=False,
positive_definite=False,
auto=False):
r"""Return `(mtype, A)` for the matrix `A`.
Parameters
----------
A : sparse array
The matrix. This is used for `A.dtype` and if `auto` is
`True`.
structurally_symmetric : bool
Set to `True` if the sparsity pattern is symmetric.
symmetric : bool
Set to `True` if the matrix is symmetric. If so, only the
upper triangular portion will be returned.
hermitian : bool
Set to `True` if the matrix is complex and Hermitian. If
so, only the upper triangular portion will be returned.
positive_definite : bool
Set to `True` if the matrix is real and symmetric positive
definite or complex and Hermitian positive definite. If
so, only the upper triangular portion will be returned.
auto : bool
If `True`, Attempt to figure out the best structure for the
matrix. Will assume not positive definite.
Returns
-------
mtype : int
Matrix type for pardiso library.
A : csr_matrix
Array packed in appropriate
:class:`scipy.sparse.csr_matrix` format. Only upper
triangular portion is non-zero if Hermitian or symmetric.
"""
if auto:
if (0 == (A - A.H).nnz):
hermitian = True
elif (0 == (A - A.T).nnz):
symmetric = True
else:
data = A.data
A.data = 0*data + 1
if (0 == (A - A.T).nnz):
structurally_symmetric == True
A.data = data
if A.dtype == float:
if positive_definite:
mtype = 2
elif symmetric:
mtype = -2
elif structurally_symmetric:
mtype = 1
else:
mtype = 11
elif A.dtype == complex:
if positive_definite:
mtype = 4
elif hermitian:
mtype = -4
elif symmetric:
mtype = 6
elif structurally_symmetric:
mtype = 3
else:
mtype = 13
else:
raise ValueError("Matrix type %s not supported." %
(str(A.dtype),))
if symmetric or hermitian or positive_definite:
A = (sp.sparse.triu(A, format='csr') +
sp.sparse.eye(*A.shape)*np.finfo(float).tiny)
return mtype, A
[docs] def solve(self, A, b, **kwargs):
r"""Return solution `x` to `A*x = b`.
Parameters
----------
A : array
Matrix. Will be converted to
:class:`scipy.sparse.csr_matrix` so provide in this format
if possible.
b : array
Right-hand side(s). If this has two-dimensions, then
`b.shape[1]` equations will be solved simultaneously.
**kwargs :
See also the parameters of :meth:`mtype` for additional
keyword arguments allowed for specifying the type of
matrix.
Returns
-------
x : array
Solution to `A*x = b`.
"""
if not sp.sparse.issparse(A):
A = sp.sparse.csr_matrix(A)
if A.shape[0] != A.shape[1]:
raise ValueError("Matrix must be square")
mtype, A = self.mtype(A=A, **kwargs)
dtype = A.dtype
x = np.empty(b.shape, dtype=dtype)
nrhs = 1
if 1 < len(b.shape):
nrhs = b.shape[1]
error = np.zeros(2, dtype='i') # Need 'i' to get an int here
# rather than long. See numpy
# ticket #1246.
args = dict(n=int(A.shape[0]),
ia=A.indptr.astype(int),
ja=A.indices.astype(int),
a=A.data,
b=b,
x=x,
nnz=int(A.nnz),
mtype=int(mtype),
nrhs=int(nrhs),
error=error,
num_procs=int(self.threads),
solver=int(self.solver),
msglvl=int(self.verbose),
)
sp.weave.inline(
code=code % dict(hash=hash(support_code)),
support_code=support_code,
arg_names=args.keys(),
local_dict=args,
compiler='gcc',
libraries=['pardiso_GNU_MACOSX32'],
verbose=2,
)
if 0 == error[0]:
msg = None
elif 1 == error[1]:
if -10 == error[0]:
msg = "Error %i: No license file found"
elif -11 == error[0]:
msg = "Error %i: License is expired: %i"
elif -12 == error[0]:
msg = "Error %i: Wrong username or hostname"
else:
msg = "Error %i: Unknown initialization error"
elif 2 == error[1]:
msg = "Error %i: during symbolic factorization"
elif 3 == error[1]:
msg = "Error %i: during numerical factorization"
elif 4 == error[1]:
msg = "Error %i: during solution"
if msg:
err = PardisoError(msg % (error[0],))
err.error = error
raise err
return x
support_code = r"""
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/* Change this if your Fortran compiler does not append underscores. */
/* e.g. the AIX compiler: #define F77_FUNC(func) func */
#ifdef AIX
#define F77_FUNC(func) func
#else
#define F77_FUNC(func) func ## _
#endif
#ifdef __cplusplus
extern "C" {
#endif
/* PARDISO prototype. */
extern int F77_FUNC(pardisoinit)
(void *, int *, int *, int *, double *, int *);
extern int F77_FUNC(pardiso)
(void *, int *, int *, int *, int *, int *,
double *, int *, int *, int *, int *, int *,
int *, double *, double *, int *, double *);
#ifdef __cplusplus
}
#endif
"""
code = r"""
/* hash(support_code) = %(hash)i to force recompile */
/* Internal solver memory pointer pt, */
/* 32-bit: int pt[64]; 64-bit: long int pt[64] */
/* or void *pt[64] should be OK on both architectures */
void *pt[64];
/* Pardiso control parameters. */
int iparm[64];
double dparm[64];
int maxfct, mnum, phase;
/* Auxiliary variables. */
char *var;
int i;
double ddum; /* Double dummy */
int idum; /* Integer dummy. */
/* -------------------------------------------------------------------- */
/* .. Setup Pardiso control parameters. */
/* -------------------------------------------------------------------- */
error[0] = 0;
F77_FUNC(pardisoinit) (pt, &mtype, &solver, iparm, dparm, &error[0]);
if (error[0] != 0) {
error[1] = 1; /* Error checking license. */
throw 1;
}
iparm[2] = num_procs;
maxfct = 1; /* Maximum number of numerical factorizations. */
mnum = 1; /* Which factorization to use. */
error[0] = 0; /* Initialize error flag */
/* -------------------------------------------------------------------- */
/* .. Convert matrix from 0-based C-notation to Fortran 1-based */
/* notation. */
/* -------------------------------------------------------------------- */
for (i = 0; i < n+1; i++) {
ia[i] += 1;
}
for (i = 0; i < nnz; i++) {
ja[i] += 1;
}
/* -------------------------------------------------------------------- */
/* .. Reordering and Symbolic Factorization. This step also allocates */
/* all memory that is necessary for the factorization. */
/* -------------------------------------------------------------------- */
phase = 11;
F77_FUNC(pardiso) (pt, &maxfct, &mnum, &mtype, &phase,
&n, a, ia, ja, &idum, &nrhs,
iparm, &msglvl, &ddum, &ddum, &error[0], dparm);
if (error[0] != 0) {
error[1] = 2; /* Error during symbolic factorization */
throw 2;
}
if (0 < msglvl) {
printf("\nReordering completed ... ");
printf("\nNumber of nonzeros in factors = %%d", iparm[17]);
printf("\nNumber of factorization MFLOPS = %%d", iparm[18]);
}
/* -------------------------------------------------------------------- */
/* .. Numerical factorization. */
/* -------------------------------------------------------------------- */
phase = 22;
F77_FUNC(pardiso) (pt, &maxfct, &mnum, &mtype, &phase,
&n, a, ia, ja, &idum, &nrhs,
iparm, &msglvl, &ddum, &ddum, &error[0], dparm);
if (error[0] != 0) {
error[1] = 3; /* ERROR during numerical factorization */
throw 3;
}
if (0 < msglvl) {
printf("\nFactorization completed ...\n ");
}
/* -------------------------------------------------------------------- */
/* .. Back substitution and iterative refinement. */
/* -------------------------------------------------------------------- */
phase = 33;
iparm[7] = 1; /* Max numbers of iterative refinement steps. */
F77_FUNC(pardiso) (pt, &maxfct, &mnum, &mtype, &phase,
&n, a, ia, ja, &idum, &nrhs,
iparm, &msglvl, b, x, &error[0], dparm);
if (error[0] != 0) {
error[1] = 4; /* ERROR during solution */
throw 4;
}
if (0 < msglvl) {
printf("\nSolve completed ...\n");
}
/* -------------------------------------------------------------------- */
/* .. Convert matrix back to 0-based C-notation. */
/* -------------------------------------------------------------------- */
for (i = 0; i < n+1; i++) {
ia[i] -= 1;
}
for (i = 0; i < nnz; i++) {
ja[i] -= 1;
}
/* -------------------------------------------------------------------- */
/* .. Termination and release of memory. */
/* -------------------------------------------------------------------- */
phase = -1; /* Release internal memory. */
F77_FUNC(pardiso) (pt, &maxfct, &mnum, &mtype, &phase,
&n, &ddum, ia, ja, &idum, &nrhs,
iparm, &msglvl, &ddum, &ddum, &error[0], dparm);
"""