Hide keyboard shortcuts

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 

2 

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 

10 

11from . import _superlu 

12 

13noScikit = False 

14try: 

15 import scikits.umfpack as umfpack 

16except ImportError: 

17 noScikit = True 

18 

19useUmfpack = not noScikit 

20 

21__all__ = ['use_solver', 'spsolve', 'splu', 'spilu', 'factorized', 

22 'MatrixRankWarning', 'spsolve_triangular'] 

23 

24 

25class MatrixRankWarning(UserWarning): 

26 pass 

27 

28 

29def use_solver(**kwargs): 

30 """ 

31 Select default sparse direct solver to be used. 

32 

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 

42 

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. 

49 

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. 

53 

54 """ 

55 if 'useUmfpack' in kwargs: 

56 globals()['useUmfpack'] = kwargs['useUmfpack'] 

57 if useUmfpack and 'assumeSortedIndices' in kwargs: 

58 umfpack.configure(assumeSortedIndices=kwargs['assumeSortedIndices']) 

59 

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 } 

68 

69 f_type = np.sctypeDict[A.dtype.name] 

70 i_type = np.sctypeDict[A.indices.dtype.name] 

71 

72 try: 

73 family = _families[(f_type, i_type)] 

74 

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) 

80 

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) 

88 

89 return family, A_new 

90 

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. 

93 

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') 

104 

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. 

112 

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]) 

119 

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. 

127 

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 """ 

138 

139 if is_pydata_spmatrix(A): 

140 A = A.to_scipy_sparse().tocsc() 

141 

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) 

146 

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)) 

152 

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) 

161 

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),)) 

166 

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

168 raise ValueError("matrix - rhs dimension mismatch (%s - %s)" 

169 % (A.shape, b.shape[0])) 

170 

171 use_umfpack = use_umfpack and useUmfpack 

172 

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() 

179 

180 if noScikit: 

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

182 

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") 

186 

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 

195 

196 if not b_is_sparse: 

197 if isspmatrix_csc(A): 

198 flag = 1 # CSC format 

199 else: 

200 flag = 0 # CSR format 

201 

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) 

213 

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) 

218 

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) 

237 

238 if is_pydata_spmatrix(b): 

239 x = b.__class__(x) 

240 

241 return x 

242 

243 

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. 

248 

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') 

256 

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 

261 

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. 

277 

278 Returns 

279 ------- 

280 invA : scipy.sparse.linalg.SuperLU 

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

282 

283 See also 

284 -------- 

285 spilu : incomplete LU decomposition 

286 

287 Notes 

288 ----- 

289 This function uses the SuperLU library. 

290 

291 References 

292 ---------- 

293 .. [1] SuperLU http://crd.lbl.gov/~xiaoye/SuperLU/ 

294 

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 """ 

309 

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 

315 

316 if not isspmatrix_csc(A): 

317 A = csc_matrix(A) 

318 warn('splu requires CSC matrix format', SparseEfficiencyWarning) 

319 

320 # sum duplicates for non-canonical format 

321 A.sum_duplicates() 

322 A = A.asfptype() # upcast to a floating point format 

323 

324 M, N = A.shape 

325 if (M != N): 

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

327 

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) 

332 

333 # Ensure that no column permutations are applied 

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

335 _options["SymmetricMode"] = True 

336 

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) 

340 

341 

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. 

346 

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

348 

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``) 

362 

363 See SuperLU documentation for details. 

364 

365 Remaining other options 

366 Same as for `splu` 

367 

368 Returns 

369 ------- 

370 invA_approx : scipy.sparse.linalg.SuperLU 

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

372 

373 See also 

374 -------- 

375 splu : complete LU decomposition 

376 

377 Notes 

378 ----- 

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

380 increase `fill_factor` AND decrease `drop_tol`. 

381 

382 This function uses the SuperLU library. 

383 

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 """ 

398 

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 

404 

405 if not isspmatrix_csc(A): 

406 A = csc_matrix(A) 

407 warn('splu requires CSC matrix format', SparseEfficiencyWarning) 

408 

409 # sum duplicates for non-canonical format 

410 A.sum_duplicates() 

411 A = A.asfptype() # upcast to a floating point format 

412 

413 M, N = A.shape 

414 if (M != N): 

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

416 

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) 

423 

424 # Ensure that no column permutations are applied 

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

426 _options["SymmetricMode"] = True 

427 

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) 

431 

432 

433def factorized(A): 

434 """ 

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

436 

437 Parameters 

438 ---------- 

439 A : (N, N) array_like 

440 Input. 

441 

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,). 

447 

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.]) 

458 

459 """ 

460 if is_pydata_spmatrix(A): 

461 A = A.to_scipy_sparse().tocsc() 

462 

463 if useUmfpack: 

464 if noScikit: 

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

466 

467 if not isspmatrix_csc(A): 

468 A = csc_matrix(A) 

469 warn('splu requires CSC matrix format', SparseEfficiencyWarning) 

470 

471 A = A.asfptype() # upcast to a floating point format 

472 

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") 

476 

477 umf_family, A = _get_umf_family(A) 

478 umf = umfpack.UmfpackContext(umf_family) 

479 

480 # Make LU decomposition. 

481 umf.numeric(A) 

482 

483 def solve(b): 

484 return umf.solve(umfpack.UMFPACK_A, A, b, autoTranspose=True) 

485 

486 return solve 

487 else: 

488 return splu(A).solve 

489 

490 

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. 

495 

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. 

517 

518 .. versionadded:: 1.4.0 

519 

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`. 

524 

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. 

531 

532 Notes 

533 ----- 

534 .. versionadded:: 0.19.0 

535 

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 """ 

546 

547 if is_pydata_spmatrix(A): 

548 A = A.to_scipy_sparse().tocsr() 

549 

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() 

557 

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)) 

561 

562 # sum duplicates for non-canonical format 

563 A.sum_duplicates() 

564 

565 b = np.asanyarray(b) 

566 

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)) 

575 

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) 

587 

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) 

593 

594 # Fill x iteratively. 

595 for i in row_indices: 

596 

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) 

606 

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])) 

616 

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) 

621 

622 # Compute i-th entry of x. 

623 if not unit_diagonal: 

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

625 

626 return x