Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/sparse/linalg/_dsolve/linsolve.py: 11%
209 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
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 [1]_, [2]_, [3]_, [4]_. over SuperLU. Has effect only
37 if ``scikits.umfpack`` is 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
41 installed. 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 References
55 ----------
56 .. [1] T. A. Davis, Algorithm 832: UMFPACK - an unsymmetric-pattern
57 multifrontal method with a column pre-ordering strategy, ACM
58 Trans. on Mathematical Software, 30(2), 2004, pp. 196--199.
59 https://dl.acm.org/doi/abs/10.1145/992200.992206
61 .. [2] T. A. Davis, A column pre-ordering strategy for the
62 unsymmetric-pattern multifrontal method, ACM Trans.
63 on Mathematical Software, 30(2), 2004, pp. 165--195.
64 https://dl.acm.org/doi/abs/10.1145/992200.992205
66 .. [3] T. A. Davis and I. S. Duff, A combined unifrontal/multifrontal
67 method for unsymmetric sparse matrices, ACM Trans. on
68 Mathematical Software, 25(1), 1999, pp. 1--19.
69 https://doi.org/10.1145/305658.287640
71 .. [4] T. A. Davis and I. S. Duff, An unsymmetric-pattern multifrontal
72 method for sparse LU factorization, SIAM J. Matrix Analysis and
73 Computations, 18(1), 1997, pp. 140--158.
74 https://doi.org/10.1137/S0895479894246905T.
76 Examples
77 --------
78 >>> import numpy as np
79 >>> from scipy.sparse.linalg import use_solver, spsolve
80 >>> from scipy.sparse import csc_matrix
81 >>> R = np.random.randn(5, 5)
82 >>> A = csc_matrix(R)
83 >>> b = np.random.randn(5)
84 >>> use_solver(useUmfpack=False) # enforce superLU over UMFPACK
85 >>> x = spsolve(A, b)
86 >>> np.allclose(A.dot(x), b)
87 True
88 >>> use_solver(useUmfpack=True) # reset umfPack usage to default
89 """
90 if 'useUmfpack' in kwargs:
91 globals()['useUmfpack'] = kwargs['useUmfpack']
92 if useUmfpack and 'assumeSortedIndices' in kwargs:
93 umfpack.configure(assumeSortedIndices=kwargs['assumeSortedIndices'])
95def _get_umf_family(A):
96 """Get umfpack family string given the sparse matrix dtype."""
97 _families = {
98 (np.float64, np.int32): 'di',
99 (np.complex128, np.int32): 'zi',
100 (np.float64, np.int64): 'dl',
101 (np.complex128, np.int64): 'zl'
102 }
104 f_type = np.sctypeDict[A.dtype.name]
105 i_type = np.sctypeDict[A.indices.dtype.name]
107 try:
108 family = _families[(f_type, i_type)]
110 except KeyError as e:
111 msg = 'only float64 or complex128 matrices with int32 or int64' \
112 ' indices are supported! (got: matrix: %s, indices: %s)' \
113 % (f_type, i_type)
114 raise ValueError(msg) from e
116 # See gh-8278. Considered converting only if
117 # A.shape[0]*A.shape[1] > np.iinfo(np.int32).max,
118 # but that didn't always fix the issue.
119 family = family[0] + "l"
120 A_new = copy.copy(A)
121 A_new.indptr = np.array(A.indptr, copy=False, dtype=np.int64)
122 A_new.indices = np.array(A.indices, copy=False, dtype=np.int64)
124 return family, A_new
126def spsolve(A, b, permc_spec=None, use_umfpack=True):
127 """Solve the sparse linear system Ax=b, where b may be a vector or a matrix.
129 Parameters
130 ----------
131 A : ndarray or sparse matrix
132 The square matrix A will be converted into CSC or CSR form
133 b : ndarray or sparse matrix
134 The matrix or vector representing the right hand side of the equation.
135 If a vector, b.shape must be (n,) or (n, 1).
136 permc_spec : str, optional
137 How to permute the columns of the matrix for sparsity preservation.
138 (default: 'COLAMD')
140 - ``NATURAL``: natural ordering.
141 - ``MMD_ATA``: minimum degree ordering on the structure of A^T A.
142 - ``MMD_AT_PLUS_A``: minimum degree ordering on the structure of A^T+A.
143 - ``COLAMD``: approximate minimum degree column ordering [1]_, [2]_.
145 use_umfpack : bool, optional
146 if True (default) then use UMFPACK for the solution [3]_, [4]_, [5]_,
147 [6]_ . This is only referenced if b is a vector and
148 ``scikits.umfpack`` is installed.
150 Returns
151 -------
152 x : ndarray or sparse matrix
153 the solution of the sparse linear equation.
154 If b is a vector, then x is a vector of size A.shape[1]
155 If b is a matrix, then x is a matrix of size (A.shape[1], b.shape[1])
157 Notes
158 -----
159 For solving the matrix expression AX = B, this solver assumes the resulting
160 matrix X is sparse, as is often the case for very sparse inputs. If the
161 resulting X is dense, the construction of this sparse result will be
162 relatively expensive. In that case, consider converting A to a dense
163 matrix and using scipy.linalg.solve or its variants.
165 References
166 ----------
167 .. [1] T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, Algorithm 836:
168 COLAMD, an approximate column minimum degree ordering algorithm,
169 ACM Trans. on Mathematical Software, 30(3), 2004, pp. 377--380.
170 :doi:`10.1145/1024074.1024080`
172 .. [2] T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, A column approximate
173 minimum degree ordering algorithm, ACM Trans. on Mathematical
174 Software, 30(3), 2004, pp. 353--376. :doi:`10.1145/1024074.1024079`
176 .. [3] T. A. Davis, Algorithm 832: UMFPACK - an unsymmetric-pattern
177 multifrontal method with a column pre-ordering strategy, ACM
178 Trans. on Mathematical Software, 30(2), 2004, pp. 196--199.
179 https://dl.acm.org/doi/abs/10.1145/992200.992206
181 .. [4] T. A. Davis, A column pre-ordering strategy for the
182 unsymmetric-pattern multifrontal method, ACM Trans.
183 on Mathematical Software, 30(2), 2004, pp. 165--195.
184 https://dl.acm.org/doi/abs/10.1145/992200.992205
186 .. [5] T. A. Davis and I. S. Duff, A combined unifrontal/multifrontal
187 method for unsymmetric sparse matrices, ACM Trans. on
188 Mathematical Software, 25(1), 1999, pp. 1--19.
189 https://doi.org/10.1145/305658.287640
191 .. [6] T. A. Davis and I. S. Duff, An unsymmetric-pattern multifrontal
192 method for sparse LU factorization, SIAM J. Matrix Analysis and
193 Computations, 18(1), 1997, pp. 140--158.
194 https://doi.org/10.1137/S0895479894246905T.
197 Examples
198 --------
199 >>> import numpy as np
200 >>> from scipy.sparse import csc_matrix
201 >>> from scipy.sparse.linalg import spsolve
202 >>> A = csc_matrix([[3, 2, 0], [1, -1, 0], [0, 5, 1]], dtype=float)
203 >>> B = csc_matrix([[2, 0], [-1, 0], [2, 0]], dtype=float)
204 >>> x = spsolve(A, B)
205 >>> np.allclose(A.dot(x).toarray(), B.toarray())
206 True
207 """
209 if is_pydata_spmatrix(A):
210 A = A.to_scipy_sparse().tocsc()
212 if not (isspmatrix_csc(A) or isspmatrix_csr(A)):
213 A = csc_matrix(A)
214 warn('spsolve requires A be CSC or CSR matrix format',
215 SparseEfficiencyWarning)
217 # b is a vector only if b have shape (n,) or (n, 1)
218 b_is_sparse = isspmatrix(b) or is_pydata_spmatrix(b)
219 if not b_is_sparse:
220 b = asarray(b)
221 b_is_vector = ((b.ndim == 1) or (b.ndim == 2 and b.shape[1] == 1))
223 # sum duplicates for non-canonical format
224 A.sum_duplicates()
225 A = A.asfptype() # upcast to a floating point format
226 result_dtype = np.promote_types(A.dtype, b.dtype)
227 if A.dtype != result_dtype:
228 A = A.astype(result_dtype)
229 if b.dtype != result_dtype:
230 b = b.astype(result_dtype)
232 # validate input shapes
233 M, N = A.shape
234 if (M != N):
235 raise ValueError("matrix must be square (has shape %s)" % ((M, N),))
237 if M != b.shape[0]:
238 raise ValueError("matrix - rhs dimension mismatch (%s - %s)"
239 % (A.shape, b.shape[0]))
241 use_umfpack = use_umfpack and useUmfpack
243 if b_is_vector and use_umfpack:
244 if b_is_sparse:
245 b_vec = b.toarray()
246 else:
247 b_vec = b
248 b_vec = asarray(b_vec, dtype=A.dtype).ravel()
250 if noScikit:
251 raise RuntimeError('Scikits.umfpack not installed.')
253 if A.dtype.char not in 'dD':
254 raise ValueError("convert matrix data to double, please, using"
255 " .astype(), or set linsolve.useUmfpack = False")
257 umf_family, A = _get_umf_family(A)
258 umf = umfpack.UmfpackContext(umf_family)
259 x = umf.linsolve(umfpack.UMFPACK_A, A, b_vec,
260 autoTranspose=True)
261 else:
262 if b_is_vector and b_is_sparse:
263 b = b.toarray()
264 b_is_sparse = False
266 if not b_is_sparse:
267 if isspmatrix_csc(A):
268 flag = 1 # CSC format
269 else:
270 flag = 0 # CSR format
272 options = dict(ColPerm=permc_spec)
273 x, info = _superlu.gssv(N, A.nnz, A.data, A.indices, A.indptr,
274 b, flag, options=options)
275 if info != 0:
276 warn("Matrix is exactly singular", MatrixRankWarning)
277 x.fill(np.nan)
278 if b_is_vector:
279 x = x.ravel()
280 else:
281 # b is sparse
282 Afactsolve = factorized(A)
284 if not (isspmatrix_csc(b) or is_pydata_spmatrix(b)):
285 warn('spsolve is more efficient when sparse b '
286 'is in the CSC matrix format', SparseEfficiencyWarning)
287 b = csc_matrix(b)
289 # Create a sparse output matrix by repeatedly applying
290 # the sparse factorization to solve columns of b.
291 data_segs = []
292 row_segs = []
293 col_segs = []
294 for j in range(b.shape[1]):
295 # TODO: replace this with
296 # bj = b[:, j].toarray().ravel()
297 # once 1D sparse arrays are supported.
298 # That is a slightly faster code path.
299 bj = b[:, [j]].toarray().ravel()
300 xj = Afactsolve(bj)
301 w = np.flatnonzero(xj)
302 segment_length = w.shape[0]
303 row_segs.append(w)
304 col_segs.append(np.full(segment_length, j, dtype=int))
305 data_segs.append(np.asarray(xj[w], dtype=A.dtype))
306 sparse_data = np.concatenate(data_segs)
307 sparse_row = np.concatenate(row_segs)
308 sparse_col = np.concatenate(col_segs)
309 x = A.__class__((sparse_data, (sparse_row, sparse_col)),
310 shape=b.shape, dtype=A.dtype)
312 if is_pydata_spmatrix(b):
313 x = b.__class__(x)
315 return x
318def splu(A, permc_spec=None, diag_pivot_thresh=None,
319 relax=None, panel_size=None, options=dict()):
320 """
321 Compute the LU decomposition of a sparse, square matrix.
323 Parameters
324 ----------
325 A : sparse matrix
326 Sparse matrix to factorize. Most efficient when provided in CSC
327 format. Other formats will be converted to CSC before factorization.
328 permc_spec : str, optional
329 How to permute the columns of the matrix for sparsity preservation.
330 (default: 'COLAMD')
332 - ``NATURAL``: natural ordering.
333 - ``MMD_ATA``: minimum degree ordering on the structure of A^T A.
334 - ``MMD_AT_PLUS_A``: minimum degree ordering on the structure of A^T+A.
335 - ``COLAMD``: approximate minimum degree column ordering
337 diag_pivot_thresh : float, optional
338 Threshold used for a diagonal entry to be an acceptable pivot.
339 See SuperLU user's guide for details [1]_
340 relax : int, optional
341 Expert option for customizing the degree of relaxing supernodes.
342 See SuperLU user's guide for details [1]_
343 panel_size : int, optional
344 Expert option for customizing the panel size.
345 See SuperLU user's guide for details [1]_
346 options : dict, optional
347 Dictionary containing additional expert options to SuperLU.
348 See SuperLU user guide [1]_ (section 2.4 on the 'Options' argument)
349 for more details. For example, you can specify
350 ``options=dict(Equil=False, IterRefine='SINGLE'))``
351 to turn equilibration off and perform a single iterative refinement.
353 Returns
354 -------
355 invA : scipy.sparse.linalg.SuperLU
356 Object, which has a ``solve`` method.
358 See also
359 --------
360 spilu : incomplete LU decomposition
362 Notes
363 -----
364 This function uses the SuperLU library.
366 References
367 ----------
368 .. [1] SuperLU https://portal.nersc.gov/project/sparse/superlu/
370 Examples
371 --------
372 >>> import numpy as np
373 >>> from scipy.sparse import csc_matrix
374 >>> from scipy.sparse.linalg import splu
375 >>> A = csc_matrix([[1., 0., 0.], [5., 0., 2.], [0., -1., 0.]], dtype=float)
376 >>> B = splu(A)
377 >>> x = np.array([1., 2., 3.], dtype=float)
378 >>> B.solve(x)
379 array([ 1. , -3. , -1.5])
380 >>> A.dot(B.solve(x))
381 array([ 1., 2., 3.])
382 >>> B.solve(A.dot(x))
383 array([ 1., 2., 3.])
384 """
386 if is_pydata_spmatrix(A):
387 csc_construct_func = lambda *a, cls=type(A): cls(csc_matrix(*a))
388 A = A.to_scipy_sparse().tocsc()
389 else:
390 csc_construct_func = csc_matrix
392 if not isspmatrix_csc(A):
393 A = csc_matrix(A)
394 warn('splu converted its input to CSC format', SparseEfficiencyWarning)
396 # sum duplicates for non-canonical format
397 A.sum_duplicates()
398 A = A.asfptype() # upcast to a floating point format
400 M, N = A.shape
401 if (M != N):
402 raise ValueError("can only factor square matrices") # is this true?
404 _options = dict(DiagPivotThresh=diag_pivot_thresh, ColPerm=permc_spec,
405 PanelSize=panel_size, Relax=relax)
406 if options is not None:
407 _options.update(options)
409 # Ensure that no column permutations are applied
410 if (_options["ColPerm"] == "NATURAL"):
411 _options["SymmetricMode"] = True
413 return _superlu.gstrf(N, A.nnz, A.data, A.indices, A.indptr,
414 csc_construct_func=csc_construct_func,
415 ilu=False, options=_options)
418def spilu(A, drop_tol=None, fill_factor=None, drop_rule=None, permc_spec=None,
419 diag_pivot_thresh=None, relax=None, panel_size=None, options=None):
420 """
421 Compute an incomplete LU decomposition for a sparse, square matrix.
423 The resulting object is an approximation to the inverse of `A`.
425 Parameters
426 ----------
427 A : (N, N) array_like
428 Sparse matrix to factorize. Most efficient when provided in CSC format.
429 Other formats will be converted to CSC before factorization.
430 drop_tol : float, optional
431 Drop tolerance (0 <= tol <= 1) for an incomplete LU decomposition.
432 (default: 1e-4)
433 fill_factor : float, optional
434 Specifies the fill ratio upper bound (>= 1.0) for ILU. (default: 10)
435 drop_rule : str, optional
436 Comma-separated string of drop rules to use.
437 Available rules: ``basic``, ``prows``, ``column``, ``area``,
438 ``secondary``, ``dynamic``, ``interp``. (Default: ``basic,area``)
440 See SuperLU documentation for details.
442 Remaining other options
443 Same as for `splu`
445 Returns
446 -------
447 invA_approx : scipy.sparse.linalg.SuperLU
448 Object, which has a ``solve`` method.
450 See also
451 --------
452 splu : complete LU decomposition
454 Notes
455 -----
456 To improve the better approximation to the inverse, you may need to
457 increase `fill_factor` AND decrease `drop_tol`.
459 This function uses the SuperLU library.
461 Examples
462 --------
463 >>> import numpy as np
464 >>> from scipy.sparse import csc_matrix
465 >>> from scipy.sparse.linalg import spilu
466 >>> A = csc_matrix([[1., 0., 0.], [5., 0., 2.], [0., -1., 0.]], dtype=float)
467 >>> B = spilu(A)
468 >>> x = np.array([1., 2., 3.], dtype=float)
469 >>> B.solve(x)
470 array([ 1. , -3. , -1.5])
471 >>> A.dot(B.solve(x))
472 array([ 1., 2., 3.])
473 >>> B.solve(A.dot(x))
474 array([ 1., 2., 3.])
475 """
477 if is_pydata_spmatrix(A):
478 csc_construct_func = lambda *a, cls=type(A): cls(csc_matrix(*a))
479 A = A.to_scipy_sparse().tocsc()
480 else:
481 csc_construct_func = csc_matrix
483 if not isspmatrix_csc(A):
484 A = csc_matrix(A)
485 warn('spilu converted its input to CSC format',
486 SparseEfficiencyWarning)
488 # sum duplicates for non-canonical format
489 A.sum_duplicates()
490 A = A.asfptype() # upcast to a floating point format
492 M, N = A.shape
493 if (M != N):
494 raise ValueError("can only factor square matrices") # is this true?
496 _options = dict(ILU_DropRule=drop_rule, ILU_DropTol=drop_tol,
497 ILU_FillFactor=fill_factor,
498 DiagPivotThresh=diag_pivot_thresh, ColPerm=permc_spec,
499 PanelSize=panel_size, Relax=relax)
500 if options is not None:
501 _options.update(options)
503 # Ensure that no column permutations are applied
504 if (_options["ColPerm"] == "NATURAL"):
505 _options["SymmetricMode"] = True
507 return _superlu.gstrf(N, A.nnz, A.data, A.indices, A.indptr,
508 csc_construct_func=csc_construct_func,
509 ilu=True, options=_options)
512def factorized(A):
513 """
514 Return a function for solving a sparse linear system, with A pre-factorized.
516 Parameters
517 ----------
518 A : (N, N) array_like
519 Input. A in CSC format is most efficient. A CSR format matrix will
520 be converted to CSC before factorization.
522 Returns
523 -------
524 solve : callable
525 To solve the linear system of equations given in `A`, the `solve`
526 callable should be passed an ndarray of shape (N,).
528 Examples
529 --------
530 >>> import numpy as np
531 >>> from scipy.sparse.linalg import factorized
532 >>> A = np.array([[ 3. , 2. , -1. ],
533 ... [ 2. , -2. , 4. ],
534 ... [-1. , 0.5, -1. ]])
535 >>> solve = factorized(A) # Makes LU decomposition.
536 >>> rhs1 = np.array([1, -2, 0])
537 >>> solve(rhs1) # Uses the LU factors.
538 array([ 1., -2., -2.])
540 """
541 if is_pydata_spmatrix(A):
542 A = A.to_scipy_sparse().tocsc()
544 if useUmfpack:
545 if noScikit:
546 raise RuntimeError('Scikits.umfpack not installed.')
548 if not isspmatrix_csc(A):
549 A = csc_matrix(A)
550 warn('splu converted its input to CSC format',
551 SparseEfficiencyWarning)
553 A = A.asfptype() # upcast to a floating point format
555 if A.dtype.char not in 'dD':
556 raise ValueError("convert matrix data to double, please, using"
557 " .astype(), or set linsolve.useUmfpack = False")
559 umf_family, A = _get_umf_family(A)
560 umf = umfpack.UmfpackContext(umf_family)
562 # Make LU decomposition.
563 umf.numeric(A)
565 def solve(b):
566 with np.errstate(divide="ignore", invalid="ignore"):
567 # Ignoring warnings with numpy >= 1.23.0, see gh-16523
568 result = umf.solve(umfpack.UMFPACK_A, A, b, autoTranspose=True)
570 return result
572 return solve
573 else:
574 return splu(A).solve
577def spsolve_triangular(A, b, lower=True, overwrite_A=False, overwrite_b=False,
578 unit_diagonal=False):
579 """
580 Solve the equation ``A x = b`` for `x`, assuming A is a triangular matrix.
582 Parameters
583 ----------
584 A : (M, M) sparse matrix
585 A sparse square triangular matrix. Should be in CSR format.
586 b : (M,) or (M, N) array_like
587 Right-hand side matrix in ``A x = b``
588 lower : bool, optional
589 Whether `A` is a lower or upper triangular matrix.
590 Default is lower triangular matrix.
591 overwrite_A : bool, optional
592 Allow changing `A`. The indices of `A` are going to be sorted and zero
593 entries are going to be removed.
594 Enabling gives a performance gain. Default is False.
595 overwrite_b : bool, optional
596 Allow overwriting data in `b`.
597 Enabling gives a performance gain. Default is False.
598 If `overwrite_b` is True, it should be ensured that
599 `b` has an appropriate dtype to be able to store the result.
600 unit_diagonal : bool, optional
601 If True, diagonal elements of `a` are assumed to be 1 and will not be
602 referenced.
604 .. versionadded:: 1.4.0
606 Returns
607 -------
608 x : (M,) or (M, N) ndarray
609 Solution to the system ``A x = b``. Shape of return matches shape
610 of `b`.
612 Raises
613 ------
614 LinAlgError
615 If `A` is singular or not triangular.
616 ValueError
617 If shape of `A` or shape of `b` do not match the requirements.
619 Notes
620 -----
621 .. versionadded:: 0.19.0
623 Examples
624 --------
625 >>> import numpy as np
626 >>> from scipy.sparse import csr_matrix
627 >>> from scipy.sparse.linalg import spsolve_triangular
628 >>> A = csr_matrix([[3, 0, 0], [1, -1, 0], [2, 0, 1]], dtype=float)
629 >>> B = np.array([[2, 0], [-1, 0], [2, 0]], dtype=float)
630 >>> x = spsolve_triangular(A, B)
631 >>> np.allclose(A.dot(x), B)
632 True
633 """
635 if is_pydata_spmatrix(A):
636 A = A.to_scipy_sparse().tocsr()
638 # Check the input for correct type and format.
639 if not isspmatrix_csr(A):
640 warn('CSR matrix format is required. Converting to CSR matrix.',
641 SparseEfficiencyWarning)
642 A = csr_matrix(A)
643 elif not overwrite_A:
644 A = A.copy()
646 if A.shape[0] != A.shape[1]:
647 raise ValueError(
648 'A must be a square matrix but its shape is {}.'.format(A.shape))
650 # sum duplicates for non-canonical format
651 A.sum_duplicates()
653 b = np.asanyarray(b)
655 if b.ndim not in [1, 2]:
656 raise ValueError(
657 'b must have 1 or 2 dims but its shape is {}.'.format(b.shape))
658 if A.shape[0] != b.shape[0]:
659 raise ValueError(
660 'The size of the dimensions of A must be equal to '
661 'the size of the first dimension of b but the shape of A is '
662 '{} and the shape of b is {}.'.format(A.shape, b.shape))
664 # Init x as (a copy of) b.
665 x_dtype = np.result_type(A.data, b, np.float64)
666 if overwrite_b:
667 if np.can_cast(b.dtype, x_dtype, casting='same_kind'):
668 x = b
669 else:
670 raise ValueError(
671 'Cannot overwrite b (dtype {}) with result '
672 'of type {}.'.format(b.dtype, x_dtype))
673 else:
674 x = b.astype(x_dtype, copy=True)
676 # Choose forward or backward order.
677 if lower:
678 row_indices = range(len(b))
679 else:
680 row_indices = range(len(b) - 1, -1, -1)
682 # Fill x iteratively.
683 for i in row_indices:
685 # Get indices for i-th row.
686 indptr_start = A.indptr[i]
687 indptr_stop = A.indptr[i + 1]
689 if lower:
690 A_diagonal_index_row_i = indptr_stop - 1
691 A_off_diagonal_indices_row_i = slice(indptr_start, indptr_stop - 1)
692 else:
693 A_diagonal_index_row_i = indptr_start
694 A_off_diagonal_indices_row_i = slice(indptr_start + 1, indptr_stop)
696 # Check regularity and triangularity of A.
697 if not unit_diagonal and (indptr_stop <= indptr_start
698 or A.indices[A_diagonal_index_row_i] < i):
699 raise LinAlgError(
700 'A is singular: diagonal {} is zero.'.format(i))
701 if not unit_diagonal and A.indices[A_diagonal_index_row_i] > i:
702 raise LinAlgError(
703 'A is not triangular: A[{}, {}] is nonzero.'
704 ''.format(i, A.indices[A_diagonal_index_row_i]))
706 # Incorporate off-diagonal entries.
707 A_column_indices_in_row_i = A.indices[A_off_diagonal_indices_row_i]
708 A_values_in_row_i = A.data[A_off_diagonal_indices_row_i]
709 x[i] -= np.dot(x[A_column_indices_in_row_i].T, A_values_in_row_i)
711 # Compute i-th entry of x.
712 if not unit_diagonal:
713 x[i] /= A.data[A_diagonal_index_row_i]
715 return x