Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/sparse/linalg/dsolve/linsolve.py : 12%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1from warnings import warn
3import numpy as np
4from numpy import asarray
5from scipy.sparse import (isspmatrix_csc, isspmatrix_csr, isspmatrix,
6 SparseEfficiencyWarning, csc_matrix, csr_matrix)
7from scipy.sparse.sputils import is_pydata_spmatrix
8from scipy.linalg import LinAlgError
9import copy
11from . import _superlu
13noScikit = False
14try:
15 import scikits.umfpack as umfpack
16except ImportError:
17 noScikit = True
19useUmfpack = not noScikit
21__all__ = ['use_solver', 'spsolve', 'splu', 'spilu', 'factorized',
22 'MatrixRankWarning', 'spsolve_triangular']
25class MatrixRankWarning(UserWarning):
26 pass
29def use_solver(**kwargs):
30 """
31 Select default sparse direct solver to be used.
33 Parameters
34 ----------
35 useUmfpack : bool, optional
36 Use UMFPACK over SuperLU. Has effect only if scikits.umfpack is
37 installed. Default: True
38 assumeSortedIndices : bool, optional
39 Allow UMFPACK to skip the step of sorting indices for a CSR/CSC matrix.
40 Has effect only if useUmfpack is True and scikits.umfpack is installed.
41 Default: False
43 Notes
44 -----
45 The default sparse solver is umfpack when available
46 (scikits.umfpack is installed). This can be changed by passing
47 useUmfpack = False, which then causes the always present SuperLU
48 based solver to be used.
50 Umfpack requires a CSR/CSC matrix to have sorted column/row indices. If
51 sure that the matrix fulfills this, pass ``assumeSortedIndices=True``
52 to gain some speed.
54 """
55 if 'useUmfpack' in kwargs:
56 globals()['useUmfpack'] = kwargs['useUmfpack']
57 if useUmfpack and 'assumeSortedIndices' in kwargs:
58 umfpack.configure(assumeSortedIndices=kwargs['assumeSortedIndices'])
60def _get_umf_family(A):
61 """Get umfpack family string given the sparse matrix dtype."""
62 _families = {
63 (np.float64, np.int32): 'di',
64 (np.complex128, np.int32): 'zi',
65 (np.float64, np.int64): 'dl',
66 (np.complex128, np.int64): 'zl'
67 }
69 f_type = np.sctypeDict[A.dtype.name]
70 i_type = np.sctypeDict[A.indices.dtype.name]
72 try:
73 family = _families[(f_type, i_type)]
75 except KeyError:
76 msg = 'only float64 or complex128 matrices with int32 or int64' \
77 ' indices are supported! (got: matrix: %s, indices: %s)' \
78 % (f_type, i_type)
79 raise ValueError(msg)
81 # See gh-8278. Considered converting only if
82 # A.shape[0]*A.shape[1] > np.iinfo(np.int32).max,
83 # but that didn't always fix the issue.
84 family = family[0] + "l"
85 A_new = copy.copy(A)
86 A_new.indptr = np.array(A.indptr, copy=False, dtype=np.int64)
87 A_new.indices = np.array(A.indices, copy=False, dtype=np.int64)
89 return family, A_new
91def spsolve(A, b, permc_spec=None, use_umfpack=True):
92 """Solve the sparse linear system Ax=b, where b may be a vector or a matrix.
94 Parameters
95 ----------
96 A : ndarray or sparse matrix
97 The square matrix A will be converted into CSC or CSR form
98 b : ndarray or sparse matrix
99 The matrix or vector representing the right hand side of the equation.
100 If a vector, b.shape must be (n,) or (n, 1).
101 permc_spec : str, optional
102 How to permute the columns of the matrix for sparsity preservation.
103 (default: 'COLAMD')
105 - ``NATURAL``: natural ordering.
106 - ``MMD_ATA``: minimum degree ordering on the structure of A^T A.
107 - ``MMD_AT_PLUS_A``: minimum degree ordering on the structure of A^T+A.
108 - ``COLAMD``: approximate minimum degree column ordering
109 use_umfpack : bool, optional
110 if True (default) then use umfpack for the solution. This is
111 only referenced if b is a vector and ``scikit-umfpack`` is installed.
113 Returns
114 -------
115 x : ndarray or sparse matrix
116 the solution of the sparse linear equation.
117 If b is a vector, then x is a vector of size A.shape[1]
118 If b is a matrix, then x is a matrix of size (A.shape[1], b.shape[1])
120 Notes
121 -----
122 For solving the matrix expression AX = B, this solver assumes the resulting
123 matrix X is sparse, as is often the case for very sparse inputs. If the
124 resulting X is dense, the construction of this sparse result will be
125 relatively expensive. In that case, consider converting A to a dense
126 matrix and using scipy.linalg.solve or its variants.
128 Examples
129 --------
130 >>> from scipy.sparse import csc_matrix
131 >>> from scipy.sparse.linalg import spsolve
132 >>> A = csc_matrix([[3, 2, 0], [1, -1, 0], [0, 5, 1]], dtype=float)
133 >>> B = csc_matrix([[2, 0], [-1, 0], [2, 0]], dtype=float)
134 >>> x = spsolve(A, B)
135 >>> np.allclose(A.dot(x).todense(), B.todense())
136 True
137 """
139 if is_pydata_spmatrix(A):
140 A = A.to_scipy_sparse().tocsc()
142 if not (isspmatrix_csc(A) or isspmatrix_csr(A)):
143 A = csc_matrix(A)
144 warn('spsolve requires A be CSC or CSR matrix format',
145 SparseEfficiencyWarning)
147 # b is a vector only if b have shape (n,) or (n, 1)
148 b_is_sparse = isspmatrix(b) or is_pydata_spmatrix(b)
149 if not b_is_sparse:
150 b = asarray(b)
151 b_is_vector = ((b.ndim == 1) or (b.ndim == 2 and b.shape[1] == 1))
153 # sum duplicates for non-canonical format
154 A.sum_duplicates()
155 A = A.asfptype() # upcast to a floating point format
156 result_dtype = np.promote_types(A.dtype, b.dtype)
157 if A.dtype != result_dtype:
158 A = A.astype(result_dtype)
159 if b.dtype != result_dtype:
160 b = b.astype(result_dtype)
162 # validate input shapes
163 M, N = A.shape
164 if (M != N):
165 raise ValueError("matrix must be square (has shape %s)" % ((M, N),))
167 if M != b.shape[0]:
168 raise ValueError("matrix - rhs dimension mismatch (%s - %s)"
169 % (A.shape, b.shape[0]))
171 use_umfpack = use_umfpack and useUmfpack
173 if b_is_vector and use_umfpack:
174 if b_is_sparse:
175 b_vec = b.toarray()
176 else:
177 b_vec = b
178 b_vec = asarray(b_vec, dtype=A.dtype).ravel()
180 if noScikit:
181 raise RuntimeError('Scikits.umfpack not installed.')
183 if A.dtype.char not in 'dD':
184 raise ValueError("convert matrix data to double, please, using"
185 " .astype(), or set linsolve.useUmfpack = False")
187 umf_family, A = _get_umf_family(A)
188 umf = umfpack.UmfpackContext(umf_family)
189 x = umf.linsolve(umfpack.UMFPACK_A, A, b_vec,
190 autoTranspose=True)
191 else:
192 if b_is_vector and b_is_sparse:
193 b = b.toarray()
194 b_is_sparse = False
196 if not b_is_sparse:
197 if isspmatrix_csc(A):
198 flag = 1 # CSC format
199 else:
200 flag = 0 # CSR format
202 options = dict(ColPerm=permc_spec)
203 x, info = _superlu.gssv(N, A.nnz, A.data, A.indices, A.indptr,
204 b, flag, options=options)
205 if info != 0:
206 warn("Matrix is exactly singular", MatrixRankWarning)
207 x.fill(np.nan)
208 if b_is_vector:
209 x = x.ravel()
210 else:
211 # b is sparse
212 Afactsolve = factorized(A)
214 if not (isspmatrix_csc(b) or is_pydata_spmatrix(b)):
215 warn('spsolve is more efficient when sparse b '
216 'is in the CSC matrix format', SparseEfficiencyWarning)
217 b = csc_matrix(b)
219 # Create a sparse output matrix by repeatedly applying
220 # the sparse factorization to solve columns of b.
221 data_segs = []
222 row_segs = []
223 col_segs = []
224 for j in range(b.shape[1]):
225 bj = np.asarray(b[:, j].todense()).ravel()
226 xj = Afactsolve(bj)
227 w = np.flatnonzero(xj)
228 segment_length = w.shape[0]
229 row_segs.append(w)
230 col_segs.append(np.full(segment_length, j, dtype=int))
231 data_segs.append(np.asarray(xj[w], dtype=A.dtype))
232 sparse_data = np.concatenate(data_segs)
233 sparse_row = np.concatenate(row_segs)
234 sparse_col = np.concatenate(col_segs)
235 x = A.__class__((sparse_data, (sparse_row, sparse_col)),
236 shape=b.shape, dtype=A.dtype)
238 if is_pydata_spmatrix(b):
239 x = b.__class__(x)
241 return x
244def splu(A, permc_spec=None, diag_pivot_thresh=None,
245 relax=None, panel_size=None, options=dict()):
246 """
247 Compute the LU decomposition of a sparse, square matrix.
249 Parameters
250 ----------
251 A : sparse matrix
252 Sparse matrix to factorize. Should be in CSR or CSC format.
253 permc_spec : str, optional
254 How to permute the columns of the matrix for sparsity preservation.
255 (default: 'COLAMD')
257 - ``NATURAL``: natural ordering.
258 - ``MMD_ATA``: minimum degree ordering on the structure of A^T A.
259 - ``MMD_AT_PLUS_A``: minimum degree ordering on the structure of A^T+A.
260 - ``COLAMD``: approximate minimum degree column ordering
262 diag_pivot_thresh : float, optional
263 Threshold used for a diagonal entry to be an acceptable pivot.
264 See SuperLU user's guide for details [1]_
265 relax : int, optional
266 Expert option for customizing the degree of relaxing supernodes.
267 See SuperLU user's guide for details [1]_
268 panel_size : int, optional
269 Expert option for customizing the panel size.
270 See SuperLU user's guide for details [1]_
271 options : dict, optional
272 Dictionary containing additional expert options to SuperLU.
273 See SuperLU user guide [1]_ (section 2.4 on the 'Options' argument)
274 for more details. For example, you can specify
275 ``options=dict(Equil=False, IterRefine='SINGLE'))``
276 to turn equilibration off and perform a single iterative refinement.
278 Returns
279 -------
280 invA : scipy.sparse.linalg.SuperLU
281 Object, which has a ``solve`` method.
283 See also
284 --------
285 spilu : incomplete LU decomposition
287 Notes
288 -----
289 This function uses the SuperLU library.
291 References
292 ----------
293 .. [1] SuperLU http://crd.lbl.gov/~xiaoye/SuperLU/
295 Examples
296 --------
297 >>> from scipy.sparse import csc_matrix
298 >>> from scipy.sparse.linalg import splu
299 >>> A = csc_matrix([[1., 0., 0.], [5., 0., 2.], [0., -1., 0.]], dtype=float)
300 >>> B = splu(A)
301 >>> x = np.array([1., 2., 3.], dtype=float)
302 >>> B.solve(x)
303 array([ 1. , -3. , -1.5])
304 >>> A.dot(B.solve(x))
305 array([ 1., 2., 3.])
306 >>> B.solve(A.dot(x))
307 array([ 1., 2., 3.])
308 """
310 if is_pydata_spmatrix(A):
311 csc_construct_func = lambda *a, cls=type(A): cls(csc_matrix(*a))
312 A = A.to_scipy_sparse().tocsc()
313 else:
314 csc_construct_func = csc_matrix
316 if not isspmatrix_csc(A):
317 A = csc_matrix(A)
318 warn('splu requires CSC matrix format', SparseEfficiencyWarning)
320 # sum duplicates for non-canonical format
321 A.sum_duplicates()
322 A = A.asfptype() # upcast to a floating point format
324 M, N = A.shape
325 if (M != N):
326 raise ValueError("can only factor square matrices") # is this true?
328 _options = dict(DiagPivotThresh=diag_pivot_thresh, ColPerm=permc_spec,
329 PanelSize=panel_size, Relax=relax)
330 if options is not None:
331 _options.update(options)
333 # Ensure that no column permutations are applied
334 if (_options["ColPerm"] == "NATURAL"):
335 _options["SymmetricMode"] = True
337 return _superlu.gstrf(N, A.nnz, A.data, A.indices, A.indptr,
338 csc_construct_func=csc_construct_func,
339 ilu=False, options=_options)
342def spilu(A, drop_tol=None, fill_factor=None, drop_rule=None, permc_spec=None,
343 diag_pivot_thresh=None, relax=None, panel_size=None, options=None):
344 """
345 Compute an incomplete LU decomposition for a sparse, square matrix.
347 The resulting object is an approximation to the inverse of `A`.
349 Parameters
350 ----------
351 A : (N, N) array_like
352 Sparse matrix to factorize
353 drop_tol : float, optional
354 Drop tolerance (0 <= tol <= 1) for an incomplete LU decomposition.
355 (default: 1e-4)
356 fill_factor : float, optional
357 Specifies the fill ratio upper bound (>= 1.0) for ILU. (default: 10)
358 drop_rule : str, optional
359 Comma-separated string of drop rules to use.
360 Available rules: ``basic``, ``prows``, ``column``, ``area``,
361 ``secondary``, ``dynamic``, ``interp``. (Default: ``basic,area``)
363 See SuperLU documentation for details.
365 Remaining other options
366 Same as for `splu`
368 Returns
369 -------
370 invA_approx : scipy.sparse.linalg.SuperLU
371 Object, which has a ``solve`` method.
373 See also
374 --------
375 splu : complete LU decomposition
377 Notes
378 -----
379 To improve the better approximation to the inverse, you may need to
380 increase `fill_factor` AND decrease `drop_tol`.
382 This function uses the SuperLU library.
384 Examples
385 --------
386 >>> from scipy.sparse import csc_matrix
387 >>> from scipy.sparse.linalg import spilu
388 >>> A = csc_matrix([[1., 0., 0.], [5., 0., 2.], [0., -1., 0.]], dtype=float)
389 >>> B = spilu(A)
390 >>> x = np.array([1., 2., 3.], dtype=float)
391 >>> B.solve(x)
392 array([ 1. , -3. , -1.5])
393 >>> A.dot(B.solve(x))
394 array([ 1., 2., 3.])
395 >>> B.solve(A.dot(x))
396 array([ 1., 2., 3.])
397 """
399 if is_pydata_spmatrix(A):
400 csc_construct_func = lambda *a, cls=type(A): cls(csc_matrix(*a))
401 A = A.to_scipy_sparse().tocsc()
402 else:
403 csc_construct_func = csc_matrix
405 if not isspmatrix_csc(A):
406 A = csc_matrix(A)
407 warn('splu requires CSC matrix format', SparseEfficiencyWarning)
409 # sum duplicates for non-canonical format
410 A.sum_duplicates()
411 A = A.asfptype() # upcast to a floating point format
413 M, N = A.shape
414 if (M != N):
415 raise ValueError("can only factor square matrices") # is this true?
417 _options = dict(ILU_DropRule=drop_rule, ILU_DropTol=drop_tol,
418 ILU_FillFactor=fill_factor,
419 DiagPivotThresh=diag_pivot_thresh, ColPerm=permc_spec,
420 PanelSize=panel_size, Relax=relax)
421 if options is not None:
422 _options.update(options)
424 # Ensure that no column permutations are applied
425 if (_options["ColPerm"] == "NATURAL"):
426 _options["SymmetricMode"] = True
428 return _superlu.gstrf(N, A.nnz, A.data, A.indices, A.indptr,
429 csc_construct_func=csc_construct_func,
430 ilu=True, options=_options)
433def factorized(A):
434 """
435 Return a function for solving a sparse linear system, with A pre-factorized.
437 Parameters
438 ----------
439 A : (N, N) array_like
440 Input.
442 Returns
443 -------
444 solve : callable
445 To solve the linear system of equations given in `A`, the `solve`
446 callable should be passed an ndarray of shape (N,).
448 Examples
449 --------
450 >>> from scipy.sparse.linalg import factorized
451 >>> A = np.array([[ 3. , 2. , -1. ],
452 ... [ 2. , -2. , 4. ],
453 ... [-1. , 0.5, -1. ]])
454 >>> solve = factorized(A) # Makes LU decomposition.
455 >>> rhs1 = np.array([1, -2, 0])
456 >>> solve(rhs1) # Uses the LU factors.
457 array([ 1., -2., -2.])
459 """
460 if is_pydata_spmatrix(A):
461 A = A.to_scipy_sparse().tocsc()
463 if useUmfpack:
464 if noScikit:
465 raise RuntimeError('Scikits.umfpack not installed.')
467 if not isspmatrix_csc(A):
468 A = csc_matrix(A)
469 warn('splu requires CSC matrix format', SparseEfficiencyWarning)
471 A = A.asfptype() # upcast to a floating point format
473 if A.dtype.char not in 'dD':
474 raise ValueError("convert matrix data to double, please, using"
475 " .astype(), or set linsolve.useUmfpack = False")
477 umf_family, A = _get_umf_family(A)
478 umf = umfpack.UmfpackContext(umf_family)
480 # Make LU decomposition.
481 umf.numeric(A)
483 def solve(b):
484 return umf.solve(umfpack.UMFPACK_A, A, b, autoTranspose=True)
486 return solve
487 else:
488 return splu(A).solve
491def spsolve_triangular(A, b, lower=True, overwrite_A=False, overwrite_b=False,
492 unit_diagonal=False):
493 """
494 Solve the equation `A x = b` for `x`, assuming A is a triangular matrix.
496 Parameters
497 ----------
498 A : (M, M) sparse matrix
499 A sparse square triangular matrix. Should be in CSR format.
500 b : (M,) or (M, N) array_like
501 Right-hand side matrix in `A x = b`
502 lower : bool, optional
503 Whether `A` is a lower or upper triangular matrix.
504 Default is lower triangular matrix.
505 overwrite_A : bool, optional
506 Allow changing `A`. The indices of `A` are going to be sorted and zero
507 entries are going to be removed.
508 Enabling gives a performance gain. Default is False.
509 overwrite_b : bool, optional
510 Allow overwriting data in `b`.
511 Enabling gives a performance gain. Default is False.
512 If `overwrite_b` is True, it should be ensured that
513 `b` has an appropriate dtype to be able to store the result.
514 unit_diagonal : bool, optional
515 If True, diagonal elements of `a` are assumed to be 1 and will not be
516 referenced.
518 .. versionadded:: 1.4.0
520 Returns
521 -------
522 x : (M,) or (M, N) ndarray
523 Solution to the system `A x = b`. Shape of return matches shape of `b`.
525 Raises
526 ------
527 LinAlgError
528 If `A` is singular or not triangular.
529 ValueError
530 If shape of `A` or shape of `b` do not match the requirements.
532 Notes
533 -----
534 .. versionadded:: 0.19.0
536 Examples
537 --------
538 >>> from scipy.sparse import csr_matrix
539 >>> from scipy.sparse.linalg import spsolve_triangular
540 >>> A = csr_matrix([[3, 0, 0], [1, -1, 0], [2, 0, 1]], dtype=float)
541 >>> B = np.array([[2, 0], [-1, 0], [2, 0]], dtype=float)
542 >>> x = spsolve_triangular(A, B)
543 >>> np.allclose(A.dot(x), B)
544 True
545 """
547 if is_pydata_spmatrix(A):
548 A = A.to_scipy_sparse().tocsr()
550 # Check the input for correct type and format.
551 if not isspmatrix_csr(A):
552 warn('CSR matrix format is required. Converting to CSR matrix.',
553 SparseEfficiencyWarning)
554 A = csr_matrix(A)
555 elif not overwrite_A:
556 A = A.copy()
558 if A.shape[0] != A.shape[1]:
559 raise ValueError(
560 'A must be a square matrix but its shape is {}.'.format(A.shape))
562 # sum duplicates for non-canonical format
563 A.sum_duplicates()
565 b = np.asanyarray(b)
567 if b.ndim not in [1, 2]:
568 raise ValueError(
569 'b must have 1 or 2 dims but its shape is {}.'.format(b.shape))
570 if A.shape[0] != b.shape[0]:
571 raise ValueError(
572 'The size of the dimensions of A must be equal to '
573 'the size of the first dimension of b but the shape of A is '
574 '{} and the shape of b is {}.'.format(A.shape, b.shape))
576 # Init x as (a copy of) b.
577 x_dtype = np.result_type(A.data, b, np.float)
578 if overwrite_b:
579 if np.can_cast(b.dtype, x_dtype, casting='same_kind'):
580 x = b
581 else:
582 raise ValueError(
583 'Cannot overwrite b (dtype {}) with result '
584 'of type {}.'.format(b.dtype, x_dtype))
585 else:
586 x = b.astype(x_dtype, copy=True)
588 # Choose forward or backward order.
589 if lower:
590 row_indices = range(len(b))
591 else:
592 row_indices = range(len(b) - 1, -1, -1)
594 # Fill x iteratively.
595 for i in row_indices:
597 # Get indices for i-th row.
598 indptr_start = A.indptr[i]
599 indptr_stop = A.indptr[i + 1]
600 if lower:
601 A_diagonal_index_row_i = indptr_stop - 1
602 A_off_diagonal_indices_row_i = slice(indptr_start, indptr_stop - 1)
603 else:
604 A_diagonal_index_row_i = indptr_start
605 A_off_diagonal_indices_row_i = slice(indptr_start + 1, indptr_stop)
607 # Check regularity and triangularity of A.
608 if not unit_diagonal and (indptr_stop <= indptr_start
609 or A.indices[A_diagonal_index_row_i] < i):
610 raise LinAlgError(
611 'A is singular: diagonal {} is zero.'.format(i))
612 if A.indices[A_diagonal_index_row_i] > i:
613 raise LinAlgError(
614 'A is not triangular: A[{}, {}] is nonzero.'
615 ''.format(i, A.indices[A_diagonal_index_row_i]))
617 # Incorporate off-diagonal entries.
618 A_column_indices_in_row_i = A.indices[A_off_diagonal_indices_row_i]
619 A_values_in_row_i = A.data[A_off_diagonal_indices_row_i]
620 x[i] -= np.dot(x[A_column_indices_in_row_i].T, A_values_in_row_i)
622 # Compute i-th entry of x.
623 if not unit_diagonal:
624 x[i] /= A.data[A_diagonal_index_row_i]
626 return x