Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/scipy/sparse/linalg/

1from warnings import warn

3import numpy as np

4from numpy import asarray

5from scipy.sparse import (issparse,

6 SparseEfficiencyWarning, csc_matrix, csr_matrix)

7from scipy.sparse._sputils import is_pydata_spmatrix, convert_pydata_sparse_to_scipy

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 }

103

104 # A.dtype.name can only be "float64" or

105 # "complex128" in control flow

106 f_type = getattr(np, A.dtype.name)

107 # control flow may allow for more index

108 # types to get through here

109 i_type = getattr(np, A.indices.dtype.name)

110

111 try:

112 family = _families[(f_type, i_type)]

113

114 except KeyError as e:

115 msg = ('only float64 or complex128 matrices with int32 or int64 '

116 f'indices are supported! (got: matrix: {f_type}, indices: {i_type})')

117 raise ValueError(msg) from e

118

119 # See gh-8278. Considered converting only if

120 # A.shape[0]*A.shape[1] > np.iinfo(np.int32).max,

121 # but that didn't always fix the issue.

122 family = family[0] + "l"

123 A_new = copy.copy(A)

124 A_new.indptr = np.asarray(A.indptr, dtype=np.int64)

125 A_new.indices = np.asarray(A.indices, dtype=np.int64)

126

127 return family, A_new

128

129def _safe_downcast_indices(A):

130 # check for safe downcasting

131 max_value = np.iinfo(np.intc).max

132

133 if A.indptr[-1] > max_value: # indptr[-1] is max b/c indptr always sorted

134 raise ValueError("indptr values too large for SuperLU")

135

136 if max(*A.shape) > max_value: # only check large enough arrays

137 if np.any(A.indices > max_value):

138 raise ValueError("indices values too large for SuperLU")

139

140 indices = A.indices.astype(np.intc, copy=False)

141 indptr = A.indptr.astype(np.intc, copy=False)

142 return indices, indptr

143

144def spsolve(A, b, permc_spec=None, use_umfpack=True):

145 """Solve the sparse linear system Ax=b, where b may be a vector or a matrix.

146

147 Parameters

148 ----------

149 A : ndarray or sparse matrix

150 The square matrix A will be converted into CSC or CSR form

151 b : ndarray or sparse matrix

152 The matrix or vector representing the right hand side of the equation.

153 If a vector, b.shape must be (n,) or (n, 1).

154 permc_spec : str, optional

155 How to permute the columns of the matrix for sparsity preservation.

156 (default: 'COLAMD')

157

158 - ``NATURAL``: natural ordering.

159 - ``MMD_ATA``: minimum degree ordering on the structure of A^T A.

160 - ``MMD_AT_PLUS_A``: minimum degree ordering on the structure of A^T+A.

161 - ``COLAMD``: approximate minimum degree column ordering [1]_, [2]_.

162

163 use_umfpack : bool, optional

164 if True (default) then use UMFPACK for the solution [3]_, [4]_, [5]_,

165 [6]_ . This is only referenced if b is a vector and

166 ``scikits.umfpack`` is installed.

167

168 Returns

169 -------

170 x : ndarray or sparse matrix

171 the solution of the sparse linear equation.

172 If b is a vector, then x is a vector of size A.shape[1]

173 If b is a matrix, then x is a matrix of size (A.shape[1], b.shape[1])

174

175 Notes

176 -----

177 For solving the matrix expression AX = B, this solver assumes the resulting

178 matrix X is sparse, as is often the case for very sparse inputs. If the

179 resulting X is dense, the construction of this sparse result will be

180 relatively expensive. In that case, consider converting A to a dense

181 matrix and using scipy.linalg.solve or its variants.

182

183 References

184 ----------

185 .. [1] T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, Algorithm 836:

186 COLAMD, an approximate column minimum degree ordering algorithm,

187 ACM Trans. on Mathematical Software, 30(3), 2004, pp. 377--380.

188 :doi:`10.1145/1024074.1024080`

189

190 .. [2] T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, A column approximate

191 minimum degree ordering algorithm, ACM Trans. on Mathematical

192 Software, 30(3), 2004, pp. 353--376. :doi:`10.1145/1024074.1024079`

193

194 .. [3] T. A. Davis, Algorithm 832: UMFPACK - an unsymmetric-pattern

195 multifrontal method with a column pre-ordering strategy, ACM

196 Trans. on Mathematical Software, 30(2), 2004, pp. 196--199.

197 https://dl.acm.org/doi/abs/10.1145/992200.992206

198

199 .. [4] T. A. Davis, A column pre-ordering strategy for the

200 unsymmetric-pattern multifrontal method, ACM Trans.

201 on Mathematical Software, 30(2), 2004, pp. 165--195.

202 https://dl.acm.org/doi/abs/10.1145/992200.992205

203

204 .. [5] T. A. Davis and I. S. Duff, A combined unifrontal/multifrontal

205 method for unsymmetric sparse matrices, ACM Trans. on

206 Mathematical Software, 25(1), 1999, pp. 1--19.

207 https://doi.org/10.1145/305658.287640

208

209 .. [6] T. A. Davis and I. S. Duff, An unsymmetric-pattern multifrontal

210 method for sparse LU factorization, SIAM J. Matrix Analysis and

211 Computations, 18(1), 1997, pp. 140--158.

212 https://doi.org/10.1137/S0895479894246905T.

213

214

215 Examples

216 --------

217 >>> import numpy as np

218 >>> from scipy.sparse import csc_matrix

219 >>> from scipy.sparse.linalg import spsolve

220 >>> A = csc_matrix([[3, 2, 0], [1, -1, 0], [0, 5, 1]], dtype=float)

221 >>> B = csc_matrix([[2, 0], [-1, 0], [2, 0]], dtype=float)

222 >>> x = spsolve(A, B)

223 >>> np.allclose(A.dot(x).toarray(), B.toarray())

224 True

225 """

226 is_pydata_sparse = is_pydata_spmatrix(b)

227 pydata_sparse_cls = b.__class__ if is_pydata_sparse else None

228 A = convert_pydata_sparse_to_scipy(A)

229 b = convert_pydata_sparse_to_scipy(b)

230

231 if not (issparse(A) and A.format in ("csc", "csr")):

232 A = csc_matrix(A)

233 warn('spsolve requires A be CSC or CSR matrix format',

234 SparseEfficiencyWarning, stacklevel=2)

235

236 # b is a vector only if b have shape (n,) or (n, 1)

237 b_is_sparse = issparse(b)

238 if not b_is_sparse:

239 b = asarray(b)

240 b_is_vector = ((b.ndim == 1) or (b.ndim == 2 and b.shape[1] == 1))

241

242 # sum duplicates for non-canonical format

243 A.sum_duplicates()

244 A = A._asfptype() # upcast to a floating point format

245 result_dtype = np.promote_types(A.dtype, b.dtype)

246 if A.dtype != result_dtype:

247 A = A.astype(result_dtype)

248 if b.dtype != result_dtype:

249 b = b.astype(result_dtype)

250

251 # validate input shapes

252 M, N = A.shape

253 if (M != N):

254 raise ValueError(f"matrix must be square (has shape {(M, N)})")

255

256 if M != b.shape[0]:

257 raise ValueError(f"matrix - rhs dimension mismatch ({A.shape} - {b.shape[0]})")

258

259 use_umfpack = use_umfpack and useUmfpack

260

261 if b_is_vector and use_umfpack:

262 if b_is_sparse:

263 b_vec = b.toarray()

264 else:

265 b_vec = b

266 b_vec = asarray(b_vec, dtype=A.dtype).ravel()

267

268 if noScikit:

269 raise RuntimeError('Scikits.umfpack not installed.')

270

271 if A.dtype.char not in 'dD':

272 raise ValueError("convert matrix data to double, please, using"

273 " .astype(), or set linsolve.useUmfpack = False")

274

275 umf_family, A = _get_umf_family(A)

276 umf = umfpack.UmfpackContext(umf_family)

277 x = umf.linsolve(umfpack.UMFPACK_A, A, b_vec,

278 autoTranspose=True)

279 else:

280 if b_is_vector and b_is_sparse:

281 b = b.toarray()

282 b_is_sparse = False

283

284 if not b_is_sparse:

285 if A.format == "csc":

286 flag = 1 # CSC format

287 else:

288 flag = 0 # CSR format

289

290 indices = A.indices.astype(np.intc, copy=False)

291 indptr = A.indptr.astype(np.intc, copy=False)

292 options = dict(ColPerm=permc_spec)

293 x, info = _superlu.gssv(N, A.nnz, A.data, indices, indptr,

294 b, flag, options=options)

295 if info != 0:

296 warn("Matrix is exactly singular", MatrixRankWarning, stacklevel=2)

297 x.fill(np.nan)

298 if b_is_vector:

299 x = x.ravel()

300 else:

301 # b is sparse

302 Afactsolve = factorized(A)

303

304 if not (b.format == "csc" or is_pydata_spmatrix(b)):

305 warn('spsolve is more efficient when sparse b '

306 'is in the CSC matrix format',

307 SparseEfficiencyWarning, stacklevel=2)

308 b = csc_matrix(b)

309

310 # Create a sparse output matrix by repeatedly applying

311 # the sparse factorization to solve columns of b.

312 data_segs = []

313 row_segs = []

314 col_segs = []

315 for j in range(b.shape[1]):

316 # TODO: replace this with

317 # bj = b[:, j].toarray().ravel()

318 # once 1D sparse arrays are supported.

319 # That is a slightly faster code path.

320 bj = b[:, [j]].toarray().ravel()

321 xj = Afactsolve(bj)

322 w = np.flatnonzero(xj)

323 segment_length = w.shape[0]

324 row_segs.append(w)

325 col_segs.append(np.full(segment_length, j, dtype=int))

326 data_segs.append(np.asarray(xj[w], dtype=A.dtype))

327 sparse_data = np.concatenate(data_segs)

328 sparse_row = np.concatenate(row_segs)

329 sparse_col = np.concatenate(col_segs)

330 x = A.__class__((sparse_data, (sparse_row, sparse_col)),

331 shape=b.shape, dtype=A.dtype)

332

333 if is_pydata_sparse:

334 x = pydata_sparse_cls.from_scipy_sparse(x)

335

336 return x

337

338

339def splu(A, permc_spec=None, diag_pivot_thresh=None,

340 relax=None, panel_size=None, options=dict()):

341 """

342 Compute the LU decomposition of a sparse, square matrix.

343

344 Parameters

345 ----------

346 A : sparse matrix

347 Sparse matrix to factorize. Most efficient when provided in CSC

348 format. Other formats will be converted to CSC before factorization.

349 permc_spec : str, optional

350 How to permute the columns of the matrix for sparsity preservation.

351 (default: 'COLAMD')

352

353 - ``NATURAL``: natural ordering.

354 - ``MMD_ATA``: minimum degree ordering on the structure of A^T A.

355 - ``MMD_AT_PLUS_A``: minimum degree ordering on the structure of A^T+A.

356 - ``COLAMD``: approximate minimum degree column ordering

357

358 diag_pivot_thresh : float, optional

359 Threshold used for a diagonal entry to be an acceptable pivot.

360 See SuperLU user's guide for details [1]_

361 relax : int, optional

362 Expert option for customizing the degree of relaxing supernodes.

363 See SuperLU user's guide for details [1]_

364 panel_size : int, optional

365 Expert option for customizing the panel size.

366 See SuperLU user's guide for details [1]_

367 options : dict, optional

368 Dictionary containing additional expert options to SuperLU.

369 See SuperLU user guide [1]_ (section 2.4 on the 'Options' argument)

370 for more details. For example, you can specify

371 ``options=dict(Equil=False, IterRefine='SINGLE'))``

372 to turn equilibration off and perform a single iterative refinement.

373

374 Returns

375 -------

376 invA : scipy.sparse.linalg.SuperLU

377 Object, which has a ``solve`` method.

378

379 See also

380 --------

381 spilu : incomplete LU decomposition

382

383 Notes

384 -----

385 This function uses the SuperLU library.

386

387 References

388 ----------

389 .. [1] SuperLU https://portal.nersc.gov/project/sparse/superlu/

390

391 Examples

392 --------

393 >>> import numpy as np

394 >>> from scipy.sparse import csc_matrix

395 >>> from scipy.sparse.linalg import splu

396 >>> A = csc_matrix([[1., 0., 0.], [5., 0., 2.], [0., -1., 0.]], dtype=float)

397 >>> B = splu(A)

398 >>> x = np.array([1., 2., 3.], dtype=float)

399 >>> B.solve(x)

400 array([ 1. , -3. , -1.5])

401 >>> A.dot(B.solve(x))

402 array([ 1., 2., 3.])

403 >>> B.solve(A.dot(x))

404 array([ 1., 2., 3.])

405 """

406

407 if is_pydata_spmatrix(A):

408 def csc_construct_func(*a, cls=type(A)):

409 return cls.from_scipy_sparse(csc_matrix(*a))

410 A = A.to_scipy_sparse().tocsc()

411 else:

412 csc_construct_func = csc_matrix

413

414 if not (issparse(A) and A.format == "csc"):

415 A = csc_matrix(A)

416 warn('splu converted its input to CSC format',

417 SparseEfficiencyWarning, stacklevel=2)

418

419 # sum duplicates for non-canonical format

420 A.sum_duplicates()

421 A = A._asfptype() # upcast to a floating point format

422

423 M, N = A.shape

424 if (M != N):

425 raise ValueError("can only factor square matrices") # is this true?

426

427 indices, indptr = _safe_downcast_indices(A)

428

429 _options = dict(DiagPivotThresh=diag_pivot_thresh, ColPerm=permc_spec,

430 PanelSize=panel_size, Relax=relax)

431 if options is not None:

432 _options.update(options)

433

434 # Ensure that no column permutations are applied

435 if (_options["ColPerm"] == "NATURAL"):

436 _options["SymmetricMode"] = True

437

438 return _superlu.gstrf(N, A.nnz, A.data, indices, indptr,

439 csc_construct_func=csc_construct_func,

440 ilu=False, options=_options)

441

442

443def spilu(A, drop_tol=None, fill_factor=None, drop_rule=None, permc_spec=None,

444 diag_pivot_thresh=None, relax=None, panel_size=None, options=None):

445 """

446 Compute an incomplete LU decomposition for a sparse, square matrix.

447

448 The resulting object is an approximation to the inverse of `A`.

449

450 Parameters

451 ----------

452 A : (N, N) array_like

453 Sparse matrix to factorize. Most efficient when provided in CSC format.

454 Other formats will be converted to CSC before factorization.

455 drop_tol : float, optional

456 Drop tolerance (0 <= tol <= 1) for an incomplete LU decomposition.

457 (default: 1e-4)

458 fill_factor : float, optional

459 Specifies the fill ratio upper bound (>= 1.0) for ILU. (default: 10)

460 drop_rule : str, optional

461 Comma-separated string of drop rules to use.

462 Available rules: ``basic``, ``prows``, ``column``, ``area``,

463 ``secondary``, ``dynamic``, ``interp``. (Default: ``basic,area``)

464

465 See SuperLU documentation for details.

466

467 Remaining other options

468 Same as for `splu`

469

470 Returns

471 -------

472 invA_approx : scipy.sparse.linalg.SuperLU

473 Object, which has a ``solve`` method.

474

475 See also

476 --------

477 splu : complete LU decomposition

478

479 Notes

480 -----

481 To improve the better approximation to the inverse, you may need to

482 increase `fill_factor` AND decrease `drop_tol`.

483

484 This function uses the SuperLU library.

485

486 Examples

487 --------

488 >>> import numpy as np

489 >>> from scipy.sparse import csc_matrix

490 >>> from scipy.sparse.linalg import spilu

491 >>> A = csc_matrix([[1., 0., 0.], [5., 0., 2.], [0., -1., 0.]], dtype=float)

492 >>> B = spilu(A)

493 >>> x = np.array([1., 2., 3.], dtype=float)

494 >>> B.solve(x)

495 array([ 1. , -3. , -1.5])

496 >>> A.dot(B.solve(x))

497 array([ 1., 2., 3.])

498 >>> B.solve(A.dot(x))

499 array([ 1., 2., 3.])

500 """

501

502 if is_pydata_spmatrix(A):

503 def csc_construct_func(*a, cls=type(A)):

504 return cls.from_scipy_sparse(csc_matrix(*a))

505 A = A.to_scipy_sparse().tocsc()

506 else:

507 csc_construct_func = csc_matrix

508

509 if not (issparse(A) and A.format == "csc"):

510 A = csc_matrix(A)

511 warn('spilu converted its input to CSC format',

512 SparseEfficiencyWarning, stacklevel=2)

513

514 # sum duplicates for non-canonical format

515 A.sum_duplicates()

516 A = A._asfptype() # upcast to a floating point format

517

518 M, N = A.shape

519 if (M != N):

520 raise ValueError("can only factor square matrices") # is this true?

521

522 indices, indptr = _safe_downcast_indices(A)

523

524 _options = dict(ILU_DropRule=drop_rule, ILU_DropTol=drop_tol,

525 ILU_FillFactor=fill_factor,

526 DiagPivotThresh=diag_pivot_thresh, ColPerm=permc_spec,

527 PanelSize=panel_size, Relax=relax)

528 if options is not None:

529 _options.update(options)

530

531 # Ensure that no column permutations are applied

532 if (_options["ColPerm"] == "NATURAL"):

533 _options["SymmetricMode"] = True

534

535 return _superlu.gstrf(N, A.nnz, A.data, indices, indptr,

536 csc_construct_func=csc_construct_func,

537 ilu=True, options=_options)

538

539

540def factorized(A):

541 """

542 Return a function for solving a sparse linear system, with A pre-factorized.

543

544 Parameters

545 ----------

546 A : (N, N) array_like

547 Input. A in CSC format is most efficient. A CSR format matrix will

548 be converted to CSC before factorization.

549

550 Returns

551 -------

552 solve : callable

553 To solve the linear system of equations given in `A`, the `solve`

554 callable should be passed an ndarray of shape (N,).

555

556 Examples

557 --------

558 >>> import numpy as np

559 >>> from scipy.sparse.linalg import factorized

560 >>> from scipy.sparse import csc_matrix

561 >>> A = np.array([[ 3. , 2. , -1. ],

562 ... [ 2. , -2. , 4. ],

563 ... [-1. , 0.5, -1. ]])

564 >>> solve = factorized(csc_matrix(A)) # Makes LU decomposition.

565 >>> rhs1 = np.array([1, -2, 0])

566 >>> solve(rhs1) # Uses the LU factors.

567 array([ 1., -2., -2.])

568

569 """

570 if is_pydata_spmatrix(A):

571 A = A.to_scipy_sparse().tocsc()

572

573 if useUmfpack:

574 if noScikit:

575 raise RuntimeError('Scikits.umfpack not installed.')

576

577 if not (issparse(A) and A.format == "csc"):

578 A = csc_matrix(A)

579 warn('splu converted its input to CSC format',

580 SparseEfficiencyWarning, stacklevel=2)

581

582 A = A._asfptype() # upcast to a floating point format

583

584 if A.dtype.char not in 'dD':

585 raise ValueError("convert matrix data to double, please, using"

586 " .astype(), or set linsolve.useUmfpack = False")

587

588 umf_family, A = _get_umf_family(A)

589 umf = umfpack.UmfpackContext(umf_family)

590

591 # Make LU decomposition.

592 umf.numeric(A)

593

594 def solve(b):

595 with np.errstate(divide="ignore", invalid="ignore"):

596 # Ignoring warnings with numpy >= 1.23.0, see gh-16523

597 result = umf.solve(umfpack.UMFPACK_A, A, b, autoTranspose=True)

598

599 return result

600

601 return solve

602 else:

603 return splu(A).solve

604

605

606def spsolve_triangular(A, b, lower=True, overwrite_A=False, overwrite_b=False,

607 unit_diagonal=False):

608 """

609 Solve the equation ``A x = b`` for `x`, assuming A is a triangular matrix.

610

611 Parameters

612 ----------

613 A : (M, M) sparse matrix

614 A sparse square triangular matrix. Should be in CSR format.

615 b : (M,) or (M, N) array_like

616 Right-hand side matrix in ``A x = b``

617 lower : bool, optional

618 Whether `A` is a lower or upper triangular matrix.

619 Default is lower triangular matrix.

620 overwrite_A : bool, optional

621 Allow changing `A`. The indices of `A` are going to be sorted and zero

622 entries are going to be removed.

623 Enabling gives a performance gain. Default is False.

624 overwrite_b : bool, optional

625 Allow overwriting data in `b`.

626 Enabling gives a performance gain. Default is False.

627 If `overwrite_b` is True, it should be ensured that

628 `b` has an appropriate dtype to be able to store the result.

629 unit_diagonal : bool, optional

630 If True, diagonal elements of `a` are assumed to be 1 and will not be

631 referenced.

632

633 .. versionadded:: 1.4.0

634

635 Returns

636 -------

637 x : (M,) or (M, N) ndarray

638 Solution to the system ``A x = b``. Shape of return matches shape

639 of `b`.

640

641 Raises

642 ------

643 LinAlgError

644 If `A` is singular or not triangular.

645 ValueError

646 If shape of `A` or shape of `b` do not match the requirements.

647

648 Notes

649 -----

650 .. versionadded:: 0.19.0

651

652 Examples

653 --------

654 >>> import numpy as np

655 >>> from scipy.sparse import csr_matrix

656 >>> from scipy.sparse.linalg import spsolve_triangular

657 >>> A = csr_matrix([[3, 0, 0], [1, -1, 0], [2, 0, 1]], dtype=float)

658 >>> B = np.array([[2, 0], [-1, 0], [2, 0]], dtype=float)

659 >>> x = spsolve_triangular(A, B)

660 >>> np.allclose(A.dot(x), B)

661 True

662 """

663

664 if is_pydata_spmatrix(A):

665 A = A.to_scipy_sparse().tocsr()

666

667 # Check the input for correct type and format.

668 if not (issparse(A) and A.format == "csr"):

669 warn('CSR matrix format is required. Converting to CSR matrix.',

670 SparseEfficiencyWarning, stacklevel=2)

671 A = csr_matrix(A)

672 elif not overwrite_A:

673 A = A.copy()

674

675 if A.shape[0] != A.shape[1]:

676 raise ValueError(

677 f'A must be a square matrix but its shape is {A.shape}.')

678

679 # sum duplicates for non-canonical format

680 A.sum_duplicates()

681

682 b = np.asanyarray(b)

683

684 if b.ndim not in [1, 2]:

685 raise ValueError(

686 f'b must have 1 or 2 dims but its shape is {b.shape}.')

687 if A.shape[0] != b.shape[0]:

688 raise ValueError(

689 'The size of the dimensions of A must be equal to '

690 'the size of the first dimension of b but the shape of A is '

691 f'{A.shape} and the shape of b is {b.shape}.'

692 )

693

694 # Init x as (a copy of) b.

695 x_dtype = np.result_type(A.data, b, np.float64)

696 if overwrite_b:

697 if np.can_cast(b.dtype, x_dtype, casting='same_kind'):

698 x = b

699 else:

700 raise ValueError(

701 f'Cannot overwrite b (dtype {b.dtype}) with result '

702 f'of type {x_dtype}.'

703 )

704 else:

705 x = b.astype(x_dtype, copy=True)

706

707 # Choose forward or backward order.

708 if lower:

709 row_indices = range(len(b))

710 else:

711 row_indices = range(len(b) - 1, -1, -1)

712

713 # Fill x iteratively.

714 for i in row_indices:

715

716 # Get indices for i-th row.

717 indptr_start = A.indptr[i]

718 indptr_stop = A.indptr[i + 1]

719

720 if lower:

721 A_diagonal_index_row_i = indptr_stop - 1

722 A_off_diagonal_indices_row_i = slice(indptr_start, indptr_stop - 1)

723 else:

724 A_diagonal_index_row_i = indptr_start

725 A_off_diagonal_indices_row_i = slice(indptr_start + 1, indptr_stop)

726

727 # Check regularity and triangularity of A.

728 if not unit_diagonal and (indptr_stop <= indptr_start

729 or A.indices[A_diagonal_index_row_i] < i):

730 raise LinAlgError(

731 f'A is singular: diagonal {i} is zero.')

732 if not unit_diagonal and A.indices[A_diagonal_index_row_i] > i:

733 raise LinAlgError(

734 'A is not triangular: A[{}, {}] is nonzero.'

735 ''.format(i, A.indices[A_diagonal_index_row_i]))

736

737 # Incorporate off-diagonal entries.

738 A_column_indices_in_row_i = A.indices[A_off_diagonal_indices_row_i]

739 A_values_in_row_i = A.data[A_off_diagonal_indices_row_i]

740 x[i] -= np.dot(x[A_column_indices_in_row_i].T, A_values_in_row_i)

741

742 # Compute i-th entry of x.

743 if not unit_diagonal:

744 x[i] /= A.data[A_diagonal_index_row_i]

745

746 return x

Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/scipy/sparse/linalg/_dsolve/linsolve.py: 11%

227 statements