Coverage for /usr/lib/python3/dist-packages/sympy/matrices/matrices.py: 32%

960 statements  

« prev     ^ index     » next       coverage.py v7.9.1, created at 2025-06-14 15:55 +0200

1import mpmath as mp 

2from collections.abc import Callable 

3 

4 

5from sympy.core.add import Add 

6from sympy.core.basic import Basic 

7from sympy.core.function import diff 

8from sympy.core.expr import Expr 

9from sympy.core.kind import _NumberKind, UndefinedKind 

10from sympy.core.mul import Mul 

11from sympy.core.power import Pow 

12from sympy.core.singleton import S 

13from sympy.core.symbol import Dummy, Symbol, uniquely_named_symbol 

14from sympy.core.sympify import sympify, _sympify 

15from sympy.functions.combinatorial.factorials import binomial, factorial 

16from sympy.functions.elementary.complexes import re 

17from sympy.functions.elementary.exponential import exp, log 

18from sympy.functions.elementary.miscellaneous import Max, Min, sqrt 

19from sympy.functions.special.tensor_functions import KroneckerDelta, LeviCivita 

20from sympy.polys import cancel 

21from sympy.printing import sstr 

22from sympy.printing.defaults import Printable 

23from sympy.printing.str import StrPrinter 

24from sympy.utilities.iterables import flatten, NotIterable, is_sequence, reshape 

25from sympy.utilities.misc import as_int, filldedent 

26 

27from .common import ( 

28 MatrixCommon, MatrixError, NonSquareMatrixError, NonInvertibleMatrixError, 

29 ShapeError, MatrixKind, a2idx) 

30 

31from .utilities import _iszero, _is_zero_after_expand_mul, _simplify 

32 

33from .determinant import ( 

34 _find_reasonable_pivot, _find_reasonable_pivot_naive, 

35 _adjugate, _charpoly, _cofactor, _cofactor_matrix, _per, 

36 _det, _det_bareiss, _det_berkowitz, _det_LU, _minor, _minor_submatrix) 

37 

38from .reductions import _is_echelon, _echelon_form, _rank, _rref 

39from .subspaces import _columnspace, _nullspace, _rowspace, _orthogonalize 

40 

41from .eigen import ( 

42 _eigenvals, _eigenvects, 

43 _bidiagonalize, _bidiagonal_decomposition, 

44 _is_diagonalizable, _diagonalize, 

45 _is_positive_definite, _is_positive_semidefinite, 

46 _is_negative_definite, _is_negative_semidefinite, _is_indefinite, 

47 _jordan_form, _left_eigenvects, _singular_values) 

48 

49from .decompositions import ( 

50 _rank_decomposition, _cholesky, _LDLdecomposition, 

51 _LUdecomposition, _LUdecomposition_Simple, _LUdecompositionFF, 

52 _singular_value_decomposition, _QRdecomposition, _upper_hessenberg_decomposition) 

53 

54from .graph import ( 

55 _connected_components, _connected_components_decomposition, 

56 _strongly_connected_components, _strongly_connected_components_decomposition) 

57 

58from .solvers import ( 

59 _diagonal_solve, _lower_triangular_solve, _upper_triangular_solve, 

60 _cholesky_solve, _LDLsolve, _LUsolve, _QRsolve, _gauss_jordan_solve, 

61 _pinv_solve, _solve, _solve_least_squares) 

62 

63from .inverse import ( 

64 _pinv, _inv_mod, _inv_ADJ, _inv_GE, _inv_LU, _inv_CH, _inv_LDL, _inv_QR, 

65 _inv, _inv_block) 

66 

67 

68class DeferredVector(Symbol, NotIterable): 

69 """A vector whose components are deferred (e.g. for use with lambdify). 

70 

71 Examples 

72 ======== 

73 

74 >>> from sympy import DeferredVector, lambdify 

75 >>> X = DeferredVector( 'X' ) 

76 >>> X 

77 X 

78 >>> expr = (X[0] + 2, X[2] + 3) 

79 >>> func = lambdify( X, expr) 

80 >>> func( [1, 2, 3] ) 

81 (3, 6) 

82 """ 

83 

84 def __getitem__(self, i): 

85 if i == -0: 

86 i = 0 

87 if i < 0: 

88 raise IndexError('DeferredVector index out of range') 

89 component_name = '%s[%d]' % (self.name, i) 

90 return Symbol(component_name) 

91 

92 def __str__(self): 

93 return sstr(self) 

94 

95 def __repr__(self): 

96 return "DeferredVector('%s')" % self.name 

97 

98 

99class MatrixDeterminant(MatrixCommon): 

100 """Provides basic matrix determinant operations. Should not be instantiated 

101 directly. See ``determinant.py`` for their implementations.""" 

102 

103 def _eval_det_bareiss(self, iszerofunc=_is_zero_after_expand_mul): 

104 return _det_bareiss(self, iszerofunc=iszerofunc) 

105 

106 def _eval_det_berkowitz(self): 

107 return _det_berkowitz(self) 

108 

109 def _eval_det_lu(self, iszerofunc=_iszero, simpfunc=None): 

110 return _det_LU(self, iszerofunc=iszerofunc, simpfunc=simpfunc) 

111 

112 def _eval_determinant(self): # for expressions.determinant.Determinant 

113 return _det(self) 

114 

115 def adjugate(self, method="berkowitz"): 

116 return _adjugate(self, method=method) 

117 

118 def charpoly(self, x='lambda', simplify=_simplify): 

119 return _charpoly(self, x=x, simplify=simplify) 

120 

121 def cofactor(self, i, j, method="berkowitz"): 

122 return _cofactor(self, i, j, method=method) 

123 

124 def cofactor_matrix(self, method="berkowitz"): 

125 return _cofactor_matrix(self, method=method) 

126 

127 def det(self, method="bareiss", iszerofunc=None): 

128 return _det(self, method=method, iszerofunc=iszerofunc) 

129 

130 def per(self): 

131 return _per(self) 

132 

133 def minor(self, i, j, method="berkowitz"): 

134 return _minor(self, i, j, method=method) 

135 

136 def minor_submatrix(self, i, j): 

137 return _minor_submatrix(self, i, j) 

138 

139 _find_reasonable_pivot.__doc__ = _find_reasonable_pivot.__doc__ 

140 _find_reasonable_pivot_naive.__doc__ = _find_reasonable_pivot_naive.__doc__ 

141 _eval_det_bareiss.__doc__ = _det_bareiss.__doc__ 

142 _eval_det_berkowitz.__doc__ = _det_berkowitz.__doc__ 

143 _eval_det_lu.__doc__ = _det_LU.__doc__ 

144 _eval_determinant.__doc__ = _det.__doc__ 

145 adjugate.__doc__ = _adjugate.__doc__ 

146 charpoly.__doc__ = _charpoly.__doc__ 

147 cofactor.__doc__ = _cofactor.__doc__ 

148 cofactor_matrix.__doc__ = _cofactor_matrix.__doc__ 

149 det.__doc__ = _det.__doc__ 

150 per.__doc__ = _per.__doc__ 

151 minor.__doc__ = _minor.__doc__ 

152 minor_submatrix.__doc__ = _minor_submatrix.__doc__ 

153 

154 

155class MatrixReductions(MatrixDeterminant): 

156 """Provides basic matrix row/column operations. Should not be instantiated 

157 directly. See ``reductions.py`` for some of their implementations.""" 

158 

159 def echelon_form(self, iszerofunc=_iszero, simplify=False, with_pivots=False): 

160 return _echelon_form(self, iszerofunc=iszerofunc, simplify=simplify, 

161 with_pivots=with_pivots) 

162 

163 @property 

164 def is_echelon(self): 

165 return _is_echelon(self) 

166 

167 def rank(self, iszerofunc=_iszero, simplify=False): 

168 return _rank(self, iszerofunc=iszerofunc, simplify=simplify) 

169 

170 def rref(self, iszerofunc=_iszero, simplify=False, pivots=True, 

171 normalize_last=True): 

172 return _rref(self, iszerofunc=iszerofunc, simplify=simplify, 

173 pivots=pivots, normalize_last=normalize_last) 

174 

175 echelon_form.__doc__ = _echelon_form.__doc__ 

176 is_echelon.__doc__ = _is_echelon.__doc__ 

177 rank.__doc__ = _rank.__doc__ 

178 rref.__doc__ = _rref.__doc__ 

179 

180 def _normalize_op_args(self, op, col, k, col1, col2, error_str="col"): 

181 """Validate the arguments for a row/column operation. ``error_str`` 

182 can be one of "row" or "col" depending on the arguments being parsed.""" 

183 if op not in ["n->kn", "n<->m", "n->n+km"]: 

184 raise ValueError("Unknown {} operation '{}'. Valid col operations " 

185 "are 'n->kn', 'n<->m', 'n->n+km'".format(error_str, op)) 

186 

187 # define self_col according to error_str 

188 self_cols = self.cols if error_str == 'col' else self.rows 

189 

190 # normalize and validate the arguments 

191 if op == "n->kn": 

192 col = col if col is not None else col1 

193 if col is None or k is None: 

194 raise ValueError("For a {0} operation 'n->kn' you must provide the " 

195 "kwargs `{0}` and `k`".format(error_str)) 

196 if not 0 <= col < self_cols: 

197 raise ValueError("This matrix does not have a {} '{}'".format(error_str, col)) 

198 

199 elif op == "n<->m": 

200 # we need two cols to swap. It does not matter 

201 # how they were specified, so gather them together and 

202 # remove `None` 

203 cols = {col, k, col1, col2}.difference([None]) 

204 if len(cols) > 2: 

205 # maybe the user left `k` by mistake? 

206 cols = {col, col1, col2}.difference([None]) 

207 if len(cols) != 2: 

208 raise ValueError("For a {0} operation 'n<->m' you must provide the " 

209 "kwargs `{0}1` and `{0}2`".format(error_str)) 

210 col1, col2 = cols 

211 if not 0 <= col1 < self_cols: 

212 raise ValueError("This matrix does not have a {} '{}'".format(error_str, col1)) 

213 if not 0 <= col2 < self_cols: 

214 raise ValueError("This matrix does not have a {} '{}'".format(error_str, col2)) 

215 

216 elif op == "n->n+km": 

217 col = col1 if col is None else col 

218 col2 = col1 if col2 is None else col2 

219 if col is None or col2 is None or k is None: 

220 raise ValueError("For a {0} operation 'n->n+km' you must provide the " 

221 "kwargs `{0}`, `k`, and `{0}2`".format(error_str)) 

222 if col == col2: 

223 raise ValueError("For a {0} operation 'n->n+km' `{0}` and `{0}2` must " 

224 "be different.".format(error_str)) 

225 if not 0 <= col < self_cols: 

226 raise ValueError("This matrix does not have a {} '{}'".format(error_str, col)) 

227 if not 0 <= col2 < self_cols: 

228 raise ValueError("This matrix does not have a {} '{}'".format(error_str, col2)) 

229 

230 else: 

231 raise ValueError('invalid operation %s' % repr(op)) 

232 

233 return op, col, k, col1, col2 

234 

235 def _eval_col_op_multiply_col_by_const(self, col, k): 

236 def entry(i, j): 

237 if j == col: 

238 return k * self[i, j] 

239 return self[i, j] 

240 return self._new(self.rows, self.cols, entry) 

241 

242 def _eval_col_op_swap(self, col1, col2): 

243 def entry(i, j): 

244 if j == col1: 

245 return self[i, col2] 

246 elif j == col2: 

247 return self[i, col1] 

248 return self[i, j] 

249 return self._new(self.rows, self.cols, entry) 

250 

251 def _eval_col_op_add_multiple_to_other_col(self, col, k, col2): 

252 def entry(i, j): 

253 if j == col: 

254 return self[i, j] + k * self[i, col2] 

255 return self[i, j] 

256 return self._new(self.rows, self.cols, entry) 

257 

258 def _eval_row_op_swap(self, row1, row2): 

259 def entry(i, j): 

260 if i == row1: 

261 return self[row2, j] 

262 elif i == row2: 

263 return self[row1, j] 

264 return self[i, j] 

265 return self._new(self.rows, self.cols, entry) 

266 

267 def _eval_row_op_multiply_row_by_const(self, row, k): 

268 def entry(i, j): 

269 if i == row: 

270 return k * self[i, j] 

271 return self[i, j] 

272 return self._new(self.rows, self.cols, entry) 

273 

274 def _eval_row_op_add_multiple_to_other_row(self, row, k, row2): 

275 def entry(i, j): 

276 if i == row: 

277 return self[i, j] + k * self[row2, j] 

278 return self[i, j] 

279 return self._new(self.rows, self.cols, entry) 

280 

281 def elementary_col_op(self, op="n->kn", col=None, k=None, col1=None, col2=None): 

282 """Performs the elementary column operation `op`. 

283 

284 `op` may be one of 

285 

286 * ``"n->kn"`` (column n goes to k*n) 

287 * ``"n<->m"`` (swap column n and column m) 

288 * ``"n->n+km"`` (column n goes to column n + k*column m) 

289 

290 Parameters 

291 ========== 

292 

293 op : string; the elementary row operation 

294 col : the column to apply the column operation 

295 k : the multiple to apply in the column operation 

296 col1 : one column of a column swap 

297 col2 : second column of a column swap or column "m" in the column operation 

298 "n->n+km" 

299 """ 

300 

301 op, col, k, col1, col2 = self._normalize_op_args(op, col, k, col1, col2, "col") 

302 

303 # now that we've validated, we're all good to dispatch 

304 if op == "n->kn": 

305 return self._eval_col_op_multiply_col_by_const(col, k) 

306 if op == "n<->m": 

307 return self._eval_col_op_swap(col1, col2) 

308 if op == "n->n+km": 

309 return self._eval_col_op_add_multiple_to_other_col(col, k, col2) 

310 

311 def elementary_row_op(self, op="n->kn", row=None, k=None, row1=None, row2=None): 

312 """Performs the elementary row operation `op`. 

313 

314 `op` may be one of 

315 

316 * ``"n->kn"`` (row n goes to k*n) 

317 * ``"n<->m"`` (swap row n and row m) 

318 * ``"n->n+km"`` (row n goes to row n + k*row m) 

319 

320 Parameters 

321 ========== 

322 

323 op : string; the elementary row operation 

324 row : the row to apply the row operation 

325 k : the multiple to apply in the row operation 

326 row1 : one row of a row swap 

327 row2 : second row of a row swap or row "m" in the row operation 

328 "n->n+km" 

329 """ 

330 

331 op, row, k, row1, row2 = self._normalize_op_args(op, row, k, row1, row2, "row") 

332 

333 # now that we've validated, we're all good to dispatch 

334 if op == "n->kn": 

335 return self._eval_row_op_multiply_row_by_const(row, k) 

336 if op == "n<->m": 

337 return self._eval_row_op_swap(row1, row2) 

338 if op == "n->n+km": 

339 return self._eval_row_op_add_multiple_to_other_row(row, k, row2) 

340 

341 

342class MatrixSubspaces(MatrixReductions): 

343 """Provides methods relating to the fundamental subspaces of a matrix. 

344 Should not be instantiated directly. See ``subspaces.py`` for their 

345 implementations.""" 

346 

347 def columnspace(self, simplify=False): 

348 return _columnspace(self, simplify=simplify) 

349 

350 def nullspace(self, simplify=False, iszerofunc=_iszero): 

351 return _nullspace(self, simplify=simplify, iszerofunc=iszerofunc) 

352 

353 def rowspace(self, simplify=False): 

354 return _rowspace(self, simplify=simplify) 

355 

356 # This is a classmethod but is converted to such later in order to allow 

357 # assignment of __doc__ since that does not work for already wrapped 

358 # classmethods in Python 3.6. 

359 def orthogonalize(cls, *vecs, **kwargs): 

360 return _orthogonalize(cls, *vecs, **kwargs) 

361 

362 columnspace.__doc__ = _columnspace.__doc__ 

363 nullspace.__doc__ = _nullspace.__doc__ 

364 rowspace.__doc__ = _rowspace.__doc__ 

365 orthogonalize.__doc__ = _orthogonalize.__doc__ 

366 

367 orthogonalize = classmethod(orthogonalize) # type:ignore 

368 

369 

370class MatrixEigen(MatrixSubspaces): 

371 """Provides basic matrix eigenvalue/vector operations. 

372 Should not be instantiated directly. See ``eigen.py`` for their 

373 implementations.""" 

374 

375 def eigenvals(self, error_when_incomplete=True, **flags): 

376 return _eigenvals(self, error_when_incomplete=error_when_incomplete, **flags) 

377 

378 def eigenvects(self, error_when_incomplete=True, iszerofunc=_iszero, **flags): 

379 return _eigenvects(self, error_when_incomplete=error_when_incomplete, 

380 iszerofunc=iszerofunc, **flags) 

381 

382 def is_diagonalizable(self, reals_only=False, **kwargs): 

383 return _is_diagonalizable(self, reals_only=reals_only, **kwargs) 

384 

385 def diagonalize(self, reals_only=False, sort=False, normalize=False): 

386 return _diagonalize(self, reals_only=reals_only, sort=sort, 

387 normalize=normalize) 

388 

389 def bidiagonalize(self, upper=True): 

390 return _bidiagonalize(self, upper=upper) 

391 

392 def bidiagonal_decomposition(self, upper=True): 

393 return _bidiagonal_decomposition(self, upper=upper) 

394 

395 @property 

396 def is_positive_definite(self): 

397 return _is_positive_definite(self) 

398 

399 @property 

400 def is_positive_semidefinite(self): 

401 return _is_positive_semidefinite(self) 

402 

403 @property 

404 def is_negative_definite(self): 

405 return _is_negative_definite(self) 

406 

407 @property 

408 def is_negative_semidefinite(self): 

409 return _is_negative_semidefinite(self) 

410 

411 @property 

412 def is_indefinite(self): 

413 return _is_indefinite(self) 

414 

415 def jordan_form(self, calc_transform=True, **kwargs): 

416 return _jordan_form(self, calc_transform=calc_transform, **kwargs) 

417 

418 def left_eigenvects(self, **flags): 

419 return _left_eigenvects(self, **flags) 

420 

421 def singular_values(self): 

422 return _singular_values(self) 

423 

424 eigenvals.__doc__ = _eigenvals.__doc__ 

425 eigenvects.__doc__ = _eigenvects.__doc__ 

426 is_diagonalizable.__doc__ = _is_diagonalizable.__doc__ 

427 diagonalize.__doc__ = _diagonalize.__doc__ 

428 is_positive_definite.__doc__ = _is_positive_definite.__doc__ 

429 is_positive_semidefinite.__doc__ = _is_positive_semidefinite.__doc__ 

430 is_negative_definite.__doc__ = _is_negative_definite.__doc__ 

431 is_negative_semidefinite.__doc__ = _is_negative_semidefinite.__doc__ 

432 is_indefinite.__doc__ = _is_indefinite.__doc__ 

433 jordan_form.__doc__ = _jordan_form.__doc__ 

434 left_eigenvects.__doc__ = _left_eigenvects.__doc__ 

435 singular_values.__doc__ = _singular_values.__doc__ 

436 bidiagonalize.__doc__ = _bidiagonalize.__doc__ 

437 bidiagonal_decomposition.__doc__ = _bidiagonal_decomposition.__doc__ 

438 

439 

440class MatrixCalculus(MatrixCommon): 

441 """Provides calculus-related matrix operations.""" 

442 

443 def diff(self, *args, **kwargs): 

444 """Calculate the derivative of each element in the matrix. 

445 ``args`` will be passed to the ``integrate`` function. 

446 

447 Examples 

448 ======== 

449 

450 >>> from sympy import Matrix 

451 >>> from sympy.abc import x, y 

452 >>> M = Matrix([[x, y], [1, 0]]) 

453 >>> M.diff(x) 

454 Matrix([ 

455 [1, 0], 

456 [0, 0]]) 

457 

458 See Also 

459 ======== 

460 

461 integrate 

462 limit 

463 """ 

464 # XXX this should be handled here rather than in Derivative 

465 from sympy.tensor.array.array_derivatives import ArrayDerivative 

466 kwargs.setdefault('evaluate', True) 

467 deriv = ArrayDerivative(self, *args, evaluate=True) 

468 if not isinstance(self, Basic): 

469 return deriv.as_mutable() 

470 else: 

471 return deriv 

472 

473 def _eval_derivative(self, arg): 

474 return self.applyfunc(lambda x: x.diff(arg)) 

475 

476 def integrate(self, *args, **kwargs): 

477 """Integrate each element of the matrix. ``args`` will 

478 be passed to the ``integrate`` function. 

479 

480 Examples 

481 ======== 

482 

483 >>> from sympy import Matrix 

484 >>> from sympy.abc import x, y 

485 >>> M = Matrix([[x, y], [1, 0]]) 

486 >>> M.integrate((x, )) 

487 Matrix([ 

488 [x**2/2, x*y], 

489 [ x, 0]]) 

490 >>> M.integrate((x, 0, 2)) 

491 Matrix([ 

492 [2, 2*y], 

493 [2, 0]]) 

494 

495 See Also 

496 ======== 

497 

498 limit 

499 diff 

500 """ 

501 return self.applyfunc(lambda x: x.integrate(*args, **kwargs)) 

502 

503 def jacobian(self, X): 

504 """Calculates the Jacobian matrix (derivative of a vector-valued function). 

505 

506 Parameters 

507 ========== 

508 

509 ``self`` : vector of expressions representing functions f_i(x_1, ..., x_n). 

510 X : set of x_i's in order, it can be a list or a Matrix 

511 

512 Both ``self`` and X can be a row or a column matrix in any order 

513 (i.e., jacobian() should always work). 

514 

515 Examples 

516 ======== 

517 

518 >>> from sympy import sin, cos, Matrix 

519 >>> from sympy.abc import rho, phi 

520 >>> X = Matrix([rho*cos(phi), rho*sin(phi), rho**2]) 

521 >>> Y = Matrix([rho, phi]) 

522 >>> X.jacobian(Y) 

523 Matrix([ 

524 [cos(phi), -rho*sin(phi)], 

525 [sin(phi), rho*cos(phi)], 

526 [ 2*rho, 0]]) 

527 >>> X = Matrix([rho*cos(phi), rho*sin(phi)]) 

528 >>> X.jacobian(Y) 

529 Matrix([ 

530 [cos(phi), -rho*sin(phi)], 

531 [sin(phi), rho*cos(phi)]]) 

532 

533 See Also 

534 ======== 

535 

536 hessian 

537 wronskian 

538 """ 

539 if not isinstance(X, MatrixBase): 

540 X = self._new(X) 

541 # Both X and ``self`` can be a row or a column matrix, so we need to make 

542 # sure all valid combinations work, but everything else fails: 

543 if self.shape[0] == 1: 

544 m = self.shape[1] 

545 elif self.shape[1] == 1: 

546 m = self.shape[0] 

547 else: 

548 raise TypeError("``self`` must be a row or a column matrix") 

549 if X.shape[0] == 1: 

550 n = X.shape[1] 

551 elif X.shape[1] == 1: 

552 n = X.shape[0] 

553 else: 

554 raise TypeError("X must be a row or a column matrix") 

555 

556 # m is the number of functions and n is the number of variables 

557 # computing the Jacobian is now easy: 

558 return self._new(m, n, lambda j, i: self[j].diff(X[i])) 

559 

560 def limit(self, *args): 

561 """Calculate the limit of each element in the matrix. 

562 ``args`` will be passed to the ``limit`` function. 

563 

564 Examples 

565 ======== 

566 

567 >>> from sympy import Matrix 

568 >>> from sympy.abc import x, y 

569 >>> M = Matrix([[x, y], [1, 0]]) 

570 >>> M.limit(x, 2) 

571 Matrix([ 

572 [2, y], 

573 [1, 0]]) 

574 

575 See Also 

576 ======== 

577 

578 integrate 

579 diff 

580 """ 

581 return self.applyfunc(lambda x: x.limit(*args)) 

582 

583 

584# https://github.com/sympy/sympy/pull/12854 

585class MatrixDeprecated(MatrixCommon): 

586 """A class to house deprecated matrix methods.""" 

587 def berkowitz_charpoly(self, x=Dummy('lambda'), simplify=_simplify): 

588 return self.charpoly(x=x) 

589 

590 def berkowitz_det(self): 

591 """Computes determinant using Berkowitz method. 

592 

593 See Also 

594 ======== 

595 

596 det 

597 berkowitz 

598 """ 

599 return self.det(method='berkowitz') 

600 

601 def berkowitz_eigenvals(self, **flags): 

602 """Computes eigenvalues of a Matrix using Berkowitz method. 

603 

604 See Also 

605 ======== 

606 

607 berkowitz 

608 """ 

609 return self.eigenvals(**flags) 

610 

611 def berkowitz_minors(self): 

612 """Computes principal minors using Berkowitz method. 

613 

614 See Also 

615 ======== 

616 

617 berkowitz 

618 """ 

619 sign, minors = self.one, [] 

620 

621 for poly in self.berkowitz(): 

622 minors.append(sign * poly[-1]) 

623 sign = -sign 

624 

625 return tuple(minors) 

626 

627 def berkowitz(self): 

628 from sympy.matrices import zeros 

629 berk = ((1,),) 

630 if not self: 

631 return berk 

632 

633 if not self.is_square: 

634 raise NonSquareMatrixError() 

635 

636 A, N = self, self.rows 

637 transforms = [0] * (N - 1) 

638 

639 for n in range(N, 1, -1): 

640 T, k = zeros(n + 1, n), n - 1 

641 

642 R, C = -A[k, :k], A[:k, k] 

643 A, a = A[:k, :k], -A[k, k] 

644 

645 items = [C] 

646 

647 for i in range(0, n - 2): 

648 items.append(A * items[i]) 

649 

650 for i, B in enumerate(items): 

651 items[i] = (R * B)[0, 0] 

652 

653 items = [self.one, a] + items 

654 

655 for i in range(n): 

656 T[i:, i] = items[:n - i + 1] 

657 

658 transforms[k - 1] = T 

659 

660 polys = [self._new([self.one, -A[0, 0]])] 

661 

662 for i, T in enumerate(transforms): 

663 polys.append(T * polys[i]) 

664 

665 return berk + tuple(map(tuple, polys)) 

666 

667 def cofactorMatrix(self, method="berkowitz"): 

668 return self.cofactor_matrix(method=method) 

669 

670 def det_bareis(self): 

671 return _det_bareiss(self) 

672 

673 def det_LU_decomposition(self): 

674 """Compute matrix determinant using LU decomposition. 

675 

676 

677 Note that this method fails if the LU decomposition itself 

678 fails. In particular, if the matrix has no inverse this method 

679 will fail. 

680 

681 TODO: Implement algorithm for sparse matrices (SFF), 

682 http://www.eecis.udel.edu/~saunders/papers/sffge/it5.ps. 

683 

684 See Also 

685 ======== 

686 

687 

688 det 

689 det_bareiss 

690 berkowitz_det 

691 """ 

692 return self.det(method='lu') 

693 

694 def jordan_cell(self, eigenval, n): 

695 return self.jordan_block(size=n, eigenvalue=eigenval) 

696 

697 def jordan_cells(self, calc_transformation=True): 

698 P, J = self.jordan_form() 

699 return P, J.get_diag_blocks() 

700 

701 def minorEntry(self, i, j, method="berkowitz"): 

702 return self.minor(i, j, method=method) 

703 

704 def minorMatrix(self, i, j): 

705 return self.minor_submatrix(i, j) 

706 

707 def permuteBkwd(self, perm): 

708 """Permute the rows of the matrix with the given permutation in reverse.""" 

709 return self.permute_rows(perm, direction='backward') 

710 

711 def permuteFwd(self, perm): 

712 """Permute the rows of the matrix with the given permutation.""" 

713 return self.permute_rows(perm, direction='forward') 

714 

715 

716@Mul._kind_dispatcher.register(_NumberKind, MatrixKind) 

717def num_mat_mul(k1, k2): 

718 """ 

719 Return MatrixKind. The element kind is selected by recursive dispatching. 

720 Do not need to dispatch in reversed order because KindDispatcher 

721 searches for this automatically. 

722 """ 

723 # Deal with Mul._kind_dispatcher's commutativity 

724 # XXX: this function is called with either k1 or k2 as MatrixKind because 

725 # the Mul kind dispatcher is commutative. Maybe it shouldn't be. Need to 

726 # swap the args here because NumberKind does not have an element_kind 

727 # attribute. 

728 if not isinstance(k2, MatrixKind): 

729 k1, k2 = k2, k1 

730 elemk = Mul._kind_dispatcher(k1, k2.element_kind) 

731 return MatrixKind(elemk) 

732 

733 

734@Mul._kind_dispatcher.register(MatrixKind, MatrixKind) 

735def mat_mat_mul(k1, k2): 

736 """ 

737 Return MatrixKind. The element kind is selected by recursive dispatching. 

738 """ 

739 elemk = Mul._kind_dispatcher(k1.element_kind, k2.element_kind) 

740 return MatrixKind(elemk) 

741 

742 

743class MatrixBase(MatrixDeprecated, 

744 MatrixCalculus, 

745 MatrixEigen, 

746 MatrixCommon, 

747 Printable): 

748 """Base class for matrix objects.""" 

749 # Added just for numpy compatibility 

750 __array_priority__ = 11 

751 

752 is_Matrix = True 

753 _class_priority = 3 

754 _sympify = staticmethod(sympify) 

755 zero = S.Zero 

756 one = S.One 

757 

758 @property 

759 def kind(self) -> MatrixKind: 

760 elem_kinds = {e.kind for e in self.flat()} 

761 if len(elem_kinds) == 1: 

762 elemkind, = elem_kinds 

763 else: 

764 elemkind = UndefinedKind 

765 return MatrixKind(elemkind) 

766 

767 def flat(self): 

768 return [self[i, j] for i in range(self.rows) for j in range(self.cols)] 

769 

770 def __array__(self, dtype=object): 

771 from .dense import matrix2numpy 

772 return matrix2numpy(self, dtype=dtype) 

773 

774 def __len__(self): 

775 """Return the number of elements of ``self``. 

776 

777 Implemented mainly so bool(Matrix()) == False. 

778 """ 

779 return self.rows * self.cols 

780 

781 def _matrix_pow_by_jordan_blocks(self, num): 

782 from sympy.matrices import diag, MutableMatrix 

783 

784 def jordan_cell_power(jc, n): 

785 N = jc.shape[0] 

786 l = jc[0,0] 

787 if l.is_zero: 

788 if N == 1 and n.is_nonnegative: 

789 jc[0,0] = l**n 

790 elif not (n.is_integer and n.is_nonnegative): 

791 raise NonInvertibleMatrixError("Non-invertible matrix can only be raised to a nonnegative integer") 

792 else: 

793 for i in range(N): 

794 jc[0,i] = KroneckerDelta(i, n) 

795 else: 

796 for i in range(N): 

797 bn = binomial(n, i) 

798 if isinstance(bn, binomial): 

799 bn = bn._eval_expand_func() 

800 jc[0,i] = l**(n-i)*bn 

801 for i in range(N): 

802 for j in range(1, N-i): 

803 jc[j,i+j] = jc [j-1,i+j-1] 

804 

805 P, J = self.jordan_form() 

806 jordan_cells = J.get_diag_blocks() 

807 # Make sure jordan_cells matrices are mutable: 

808 jordan_cells = [MutableMatrix(j) for j in jordan_cells] 

809 for j in jordan_cells: 

810 jordan_cell_power(j, num) 

811 return self._new(P.multiply(diag(*jordan_cells)) 

812 .multiply(P.inv())) 

813 

814 def __str__(self): 

815 if S.Zero in self.shape: 

816 return 'Matrix(%s, %s, [])' % (self.rows, self.cols) 

817 return "Matrix(%s)" % str(self.tolist()) 

818 

819 def _format_str(self, printer=None): 

820 if not printer: 

821 printer = StrPrinter() 

822 # Handle zero dimensions: 

823 if S.Zero in self.shape: 

824 return 'Matrix(%s, %s, [])' % (self.rows, self.cols) 

825 if self.rows == 1: 

826 return "Matrix([%s])" % self.table(printer, rowsep=',\n') 

827 return "Matrix([\n%s])" % self.table(printer, rowsep=',\n') 

828 

829 @classmethod 

830 def irregular(cls, ntop, *matrices, **kwargs): 

831 """Return a matrix filled by the given matrices which 

832 are listed in order of appearance from left to right, top to 

833 bottom as they first appear in the matrix. They must fill the 

834 matrix completely. 

835 

836 Examples 

837 ======== 

838 

839 >>> from sympy import ones, Matrix 

840 >>> Matrix.irregular(3, ones(2,1), ones(3,3)*2, ones(2,2)*3, 

841 ... ones(1,1)*4, ones(2,2)*5, ones(1,2)*6, ones(1,2)*7) 

842 Matrix([ 

843 [1, 2, 2, 2, 3, 3], 

844 [1, 2, 2, 2, 3, 3], 

845 [4, 2, 2, 2, 5, 5], 

846 [6, 6, 7, 7, 5, 5]]) 

847 """ 

848 ntop = as_int(ntop) 

849 # make sure we are working with explicit matrices 

850 b = [i.as_explicit() if hasattr(i, 'as_explicit') else i 

851 for i in matrices] 

852 q = list(range(len(b))) 

853 dat = [i.rows for i in b] 

854 active = [q.pop(0) for _ in range(ntop)] 

855 cols = sum([b[i].cols for i in active]) 

856 rows = [] 

857 while any(dat): 

858 r = [] 

859 for a, j in enumerate(active): 

860 r.extend(b[j][-dat[j], :]) 

861 dat[j] -= 1 

862 if dat[j] == 0 and q: 

863 active[a] = q.pop(0) 

864 if len(r) != cols: 

865 raise ValueError(filldedent(''' 

866 Matrices provided do not appear to fill 

867 the space completely.''')) 

868 rows.append(r) 

869 return cls._new(rows) 

870 

871 @classmethod 

872 def _handle_ndarray(cls, arg): 

873 # NumPy array or matrix or some other object that implements 

874 # __array__. So let's first use this method to get a 

875 # numpy.array() and then make a Python list out of it. 

876 arr = arg.__array__() 

877 if len(arr.shape) == 2: 

878 rows, cols = arr.shape[0], arr.shape[1] 

879 flat_list = [cls._sympify(i) for i in arr.ravel()] 

880 return rows, cols, flat_list 

881 elif len(arr.shape) == 1: 

882 flat_list = [cls._sympify(i) for i in arr] 

883 return arr.shape[0], 1, flat_list 

884 else: 

885 raise NotImplementedError( 

886 "SymPy supports just 1D and 2D matrices") 

887 

888 @classmethod 

889 def _handle_creation_inputs(cls, *args, **kwargs): 

890 """Return the number of rows, cols and flat matrix elements. 

891 

892 Examples 

893 ======== 

894 

895 >>> from sympy import Matrix, I 

896 

897 Matrix can be constructed as follows: 

898 

899 * from a nested list of iterables 

900 

901 >>> Matrix( ((1, 2+I), (3, 4)) ) 

902 Matrix([ 

903 [1, 2 + I], 

904 [3, 4]]) 

905 

906 * from un-nested iterable (interpreted as a column) 

907 

908 >>> Matrix( [1, 2] ) 

909 Matrix([ 

910 [1], 

911 [2]]) 

912 

913 * from un-nested iterable with dimensions 

914 

915 >>> Matrix(1, 2, [1, 2] ) 

916 Matrix([[1, 2]]) 

917 

918 * from no arguments (a 0 x 0 matrix) 

919 

920 >>> Matrix() 

921 Matrix(0, 0, []) 

922 

923 * from a rule 

924 

925 >>> Matrix(2, 2, lambda i, j: i/(j + 1) ) 

926 Matrix([ 

927 [0, 0], 

928 [1, 1/2]]) 

929 

930 See Also 

931 ======== 

932 irregular - filling a matrix with irregular blocks 

933 """ 

934 from sympy.matrices import SparseMatrix 

935 from sympy.matrices.expressions.matexpr import MatrixSymbol 

936 from sympy.matrices.expressions.blockmatrix import BlockMatrix 

937 

938 flat_list = None 

939 

940 if len(args) == 1: 

941 # Matrix(SparseMatrix(...)) 

942 if isinstance(args[0], SparseMatrix): 

943 return args[0].rows, args[0].cols, flatten(args[0].tolist()) 

944 

945 # Matrix(Matrix(...)) 

946 elif isinstance(args[0], MatrixBase): 

947 return args[0].rows, args[0].cols, args[0].flat() 

948 

949 # Matrix(MatrixSymbol('X', 2, 2)) 

950 elif isinstance(args[0], Basic) and args[0].is_Matrix: 

951 return args[0].rows, args[0].cols, args[0].as_explicit().flat() 

952 

953 elif isinstance(args[0], mp.matrix): 

954 M = args[0] 

955 flat_list = [cls._sympify(x) for x in M] 

956 return M.rows, M.cols, flat_list 

957 

958 # Matrix(numpy.ones((2, 2))) 

959 elif hasattr(args[0], "__array__"): 

960 return cls._handle_ndarray(args[0]) 

961 

962 # Matrix([1, 2, 3]) or Matrix([[1, 2], [3, 4]]) 

963 elif is_sequence(args[0]) \ 

964 and not isinstance(args[0], DeferredVector): 

965 dat = list(args[0]) 

966 ismat = lambda i: isinstance(i, MatrixBase) and ( 

967 evaluate or 

968 isinstance(i, BlockMatrix) or 

969 isinstance(i, MatrixSymbol)) 

970 raw = lambda i: is_sequence(i) and not ismat(i) 

971 evaluate = kwargs.get('evaluate', True) 

972 

973 

974 if evaluate: 

975 

976 def make_explicit(x): 

977 """make Block and Symbol explicit""" 

978 if isinstance(x, BlockMatrix): 

979 return x.as_explicit() 

980 elif isinstance(x, MatrixSymbol) and all(_.is_Integer for _ in x.shape): 

981 return x.as_explicit() 

982 else: 

983 return x 

984 

985 def make_explicit_row(row): 

986 # Could be list or could be list of lists 

987 if isinstance(row, (list, tuple)): 

988 return [make_explicit(x) for x in row] 

989 else: 

990 return make_explicit(row) 

991 

992 if isinstance(dat, (list, tuple)): 

993 dat = [make_explicit_row(row) for row in dat] 

994 

995 if dat in ([], [[]]): 

996 rows = cols = 0 

997 flat_list = [] 

998 elif not any(raw(i) or ismat(i) for i in dat): 

999 # a column as a list of values 

1000 flat_list = [cls._sympify(i) for i in dat] 

1001 rows = len(flat_list) 

1002 cols = 1 if rows else 0 

1003 elif evaluate and all(ismat(i) for i in dat): 

1004 # a column as a list of matrices 

1005 ncol = {i.cols for i in dat if any(i.shape)} 

1006 if ncol: 

1007 if len(ncol) != 1: 

1008 raise ValueError('mismatched dimensions') 

1009 flat_list = [_ for i in dat for r in i.tolist() for _ in r] 

1010 cols = ncol.pop() 

1011 rows = len(flat_list)//cols 

1012 else: 

1013 rows = cols = 0 

1014 flat_list = [] 

1015 elif evaluate and any(ismat(i) for i in dat): 

1016 ncol = set() 

1017 flat_list = [] 

1018 for i in dat: 

1019 if ismat(i): 

1020 flat_list.extend( 

1021 [k for j in i.tolist() for k in j]) 

1022 if any(i.shape): 

1023 ncol.add(i.cols) 

1024 elif raw(i): 

1025 if i: 

1026 ncol.add(len(i)) 

1027 flat_list.extend([cls._sympify(ij) for ij in i]) 

1028 else: 

1029 ncol.add(1) 

1030 flat_list.append(i) 

1031 if len(ncol) > 1: 

1032 raise ValueError('mismatched dimensions') 

1033 cols = ncol.pop() 

1034 rows = len(flat_list)//cols 

1035 else: 

1036 # list of lists; each sublist is a logical row 

1037 # which might consist of many rows if the values in 

1038 # the row are matrices 

1039 flat_list = [] 

1040 ncol = set() 

1041 rows = cols = 0 

1042 for row in dat: 

1043 if not is_sequence(row) and \ 

1044 not getattr(row, 'is_Matrix', False): 

1045 raise ValueError('expecting list of lists') 

1046 

1047 if hasattr(row, '__array__'): 

1048 if 0 in row.shape: 

1049 continue 

1050 elif not row: 

1051 continue 

1052 

1053 if evaluate and all(ismat(i) for i in row): 

1054 r, c, flatT = cls._handle_creation_inputs( 

1055 [i.T for i in row]) 

1056 T = reshape(flatT, [c]) 

1057 flat = \ 

1058 [T[i][j] for j in range(c) for i in range(r)] 

1059 r, c = c, r 

1060 else: 

1061 r = 1 

1062 if getattr(row, 'is_Matrix', False): 

1063 c = 1 

1064 flat = [row] 

1065 else: 

1066 c = len(row) 

1067 flat = [cls._sympify(i) for i in row] 

1068 ncol.add(c) 

1069 if len(ncol) > 1: 

1070 raise ValueError('mismatched dimensions') 

1071 flat_list.extend(flat) 

1072 rows += r 

1073 cols = ncol.pop() if ncol else 0 

1074 

1075 elif len(args) == 3: 

1076 rows = as_int(args[0]) 

1077 cols = as_int(args[1]) 

1078 

1079 if rows < 0 or cols < 0: 

1080 raise ValueError("Cannot create a {} x {} matrix. " 

1081 "Both dimensions must be positive".format(rows, cols)) 

1082 

1083 # Matrix(2, 2, lambda i, j: i+j) 

1084 if len(args) == 3 and isinstance(args[2], Callable): 

1085 op = args[2] 

1086 flat_list = [] 

1087 for i in range(rows): 

1088 flat_list.extend( 

1089 [cls._sympify(op(cls._sympify(i), cls._sympify(j))) 

1090 for j in range(cols)]) 

1091 

1092 # Matrix(2, 2, [1, 2, 3, 4]) 

1093 elif len(args) == 3 and is_sequence(args[2]): 

1094 flat_list = args[2] 

1095 if len(flat_list) != rows * cols: 

1096 raise ValueError( 

1097 'List length should be equal to rows*columns') 

1098 flat_list = [cls._sympify(i) for i in flat_list] 

1099 

1100 

1101 # Matrix() 

1102 elif len(args) == 0: 

1103 # Empty Matrix 

1104 rows = cols = 0 

1105 flat_list = [] 

1106 

1107 if flat_list is None: 

1108 raise TypeError(filldedent(''' 

1109 Data type not understood; expecting list of lists 

1110 or lists of values.''')) 

1111 

1112 return rows, cols, flat_list 

1113 

1114 def _setitem(self, key, value): 

1115 """Helper to set value at location given by key. 

1116 

1117 Examples 

1118 ======== 

1119 

1120 >>> from sympy import Matrix, I, zeros, ones 

1121 >>> m = Matrix(((1, 2+I), (3, 4))) 

1122 >>> m 

1123 Matrix([ 

1124 [1, 2 + I], 

1125 [3, 4]]) 

1126 >>> m[1, 0] = 9 

1127 >>> m 

1128 Matrix([ 

1129 [1, 2 + I], 

1130 [9, 4]]) 

1131 >>> m[1, 0] = [[0, 1]] 

1132 

1133 To replace row r you assign to position r*m where m 

1134 is the number of columns: 

1135 

1136 >>> M = zeros(4) 

1137 >>> m = M.cols 

1138 >>> M[3*m] = ones(1, m)*2; M 

1139 Matrix([ 

1140 [0, 0, 0, 0], 

1141 [0, 0, 0, 0], 

1142 [0, 0, 0, 0], 

1143 [2, 2, 2, 2]]) 

1144 

1145 And to replace column c you can assign to position c: 

1146 

1147 >>> M[2] = ones(m, 1)*4; M 

1148 Matrix([ 

1149 [0, 0, 4, 0], 

1150 [0, 0, 4, 0], 

1151 [0, 0, 4, 0], 

1152 [2, 2, 4, 2]]) 

1153 """ 

1154 from .dense import Matrix 

1155 

1156 is_slice = isinstance(key, slice) 

1157 i, j = key = self.key2ij(key) 

1158 is_mat = isinstance(value, MatrixBase) 

1159 if isinstance(i, slice) or isinstance(j, slice): 

1160 if is_mat: 

1161 self.copyin_matrix(key, value) 

1162 return 

1163 if not isinstance(value, Expr) and is_sequence(value): 

1164 self.copyin_list(key, value) 

1165 return 

1166 raise ValueError('unexpected value: %s' % value) 

1167 else: 

1168 if (not is_mat and 

1169 not isinstance(value, Basic) and is_sequence(value)): 

1170 value = Matrix(value) 

1171 is_mat = True 

1172 if is_mat: 

1173 if is_slice: 

1174 key = (slice(*divmod(i, self.cols)), 

1175 slice(*divmod(j, self.cols))) 

1176 else: 

1177 key = (slice(i, i + value.rows), 

1178 slice(j, j + value.cols)) 

1179 self.copyin_matrix(key, value) 

1180 else: 

1181 return i, j, self._sympify(value) 

1182 return 

1183 

1184 def add(self, b): 

1185 """Return self + b.""" 

1186 return self + b 

1187 

1188 def condition_number(self): 

1189 """Returns the condition number of a matrix. 

1190 

1191 This is the maximum singular value divided by the minimum singular value 

1192 

1193 Examples 

1194 ======== 

1195 

1196 >>> from sympy import Matrix, S 

1197 >>> A = Matrix([[1, 0, 0], [0, 10, 0], [0, 0, S.One/10]]) 

1198 >>> A.condition_number() 

1199 100 

1200 

1201 See Also 

1202 ======== 

1203 

1204 singular_values 

1205 """ 

1206 

1207 if not self: 

1208 return self.zero 

1209 singularvalues = self.singular_values() 

1210 return Max(*singularvalues) / Min(*singularvalues) 

1211 

1212 def copy(self): 

1213 """ 

1214 Returns the copy of a matrix. 

1215 

1216 Examples 

1217 ======== 

1218 

1219 >>> from sympy import Matrix 

1220 >>> A = Matrix(2, 2, [1, 2, 3, 4]) 

1221 >>> A.copy() 

1222 Matrix([ 

1223 [1, 2], 

1224 [3, 4]]) 

1225 

1226 """ 

1227 return self._new(self.rows, self.cols, self.flat()) 

1228 

1229 def cross(self, b): 

1230 r""" 

1231 Return the cross product of ``self`` and ``b`` relaxing the condition 

1232 of compatible dimensions: if each has 3 elements, a matrix of the 

1233 same type and shape as ``self`` will be returned. If ``b`` has the same 

1234 shape as ``self`` then common identities for the cross product (like 

1235 `a \times b = - b \times a`) will hold. 

1236 

1237 Parameters 

1238 ========== 

1239 b : 3x1 or 1x3 Matrix 

1240 

1241 See Also 

1242 ======== 

1243 

1244 dot 

1245 multiply 

1246 multiply_elementwise 

1247 """ 

1248 from sympy.matrices.expressions.matexpr import MatrixExpr 

1249 

1250 if not isinstance(b, (MatrixBase, MatrixExpr)): 

1251 raise TypeError( 

1252 "{} must be a Matrix, not {}.".format(b, type(b))) 

1253 

1254 if not (self.rows * self.cols == b.rows * b.cols == 3): 

1255 raise ShapeError("Dimensions incorrect for cross product: %s x %s" % 

1256 ((self.rows, self.cols), (b.rows, b.cols))) 

1257 else: 

1258 return self._new(self.rows, self.cols, ( 

1259 (self[1] * b[2] - self[2] * b[1]), 

1260 (self[2] * b[0] - self[0] * b[2]), 

1261 (self[0] * b[1] - self[1] * b[0]))) 

1262 

1263 @property 

1264 def D(self): 

1265 """Return Dirac conjugate (if ``self.rows == 4``). 

1266 

1267 Examples 

1268 ======== 

1269 

1270 >>> from sympy import Matrix, I, eye 

1271 >>> m = Matrix((0, 1 + I, 2, 3)) 

1272 >>> m.D 

1273 Matrix([[0, 1 - I, -2, -3]]) 

1274 >>> m = (eye(4) + I*eye(4)) 

1275 >>> m[0, 3] = 2 

1276 >>> m.D 

1277 Matrix([ 

1278 [1 - I, 0, 0, 0], 

1279 [ 0, 1 - I, 0, 0], 

1280 [ 0, 0, -1 + I, 0], 

1281 [ 2, 0, 0, -1 + I]]) 

1282 

1283 If the matrix does not have 4 rows an AttributeError will be raised 

1284 because this property is only defined for matrices with 4 rows. 

1285 

1286 >>> Matrix(eye(2)).D 

1287 Traceback (most recent call last): 

1288 ... 

1289 AttributeError: Matrix has no attribute D. 

1290 

1291 See Also 

1292 ======== 

1293 

1294 sympy.matrices.common.MatrixCommon.conjugate: By-element conjugation 

1295 sympy.matrices.common.MatrixCommon.H: Hermite conjugation 

1296 """ 

1297 from sympy.physics.matrices import mgamma 

1298 if self.rows != 4: 

1299 # In Python 3.2, properties can only return an AttributeError 

1300 # so we can't raise a ShapeError -- see commit which added the 

1301 # first line of this inline comment. Also, there is no need 

1302 # for a message since MatrixBase will raise the AttributeError 

1303 raise AttributeError 

1304 return self.H * mgamma(0) 

1305 

1306 def dot(self, b, hermitian=None, conjugate_convention=None): 

1307 """Return the dot or inner product of two vectors of equal length. 

1308 Here ``self`` must be a ``Matrix`` of size 1 x n or n x 1, and ``b`` 

1309 must be either a matrix of size 1 x n, n x 1, or a list/tuple of length n. 

1310 A scalar is returned. 

1311 

1312 By default, ``dot`` does not conjugate ``self`` or ``b``, even if there are 

1313 complex entries. Set ``hermitian=True`` (and optionally a ``conjugate_convention``) 

1314 to compute the hermitian inner product. 

1315 

1316 Possible kwargs are ``hermitian`` and ``conjugate_convention``. 

1317 

1318 If ``conjugate_convention`` is ``"left"``, ``"math"`` or ``"maths"``, 

1319 the conjugate of the first vector (``self``) is used. If ``"right"`` 

1320 or ``"physics"`` is specified, the conjugate of the second vector ``b`` is used. 

1321 

1322 Examples 

1323 ======== 

1324 

1325 >>> from sympy import Matrix 

1326 >>> M = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) 

1327 >>> v = Matrix([1, 1, 1]) 

1328 >>> M.row(0).dot(v) 

1329 6 

1330 >>> M.col(0).dot(v) 

1331 12 

1332 >>> v = [3, 2, 1] 

1333 >>> M.row(0).dot(v) 

1334 10 

1335 

1336 >>> from sympy import I 

1337 >>> q = Matrix([1*I, 1*I, 1*I]) 

1338 >>> q.dot(q, hermitian=False) 

1339 -3 

1340 

1341 >>> q.dot(q, hermitian=True) 

1342 3 

1343 

1344 >>> q1 = Matrix([1, 1, 1*I]) 

1345 >>> q.dot(q1, hermitian=True, conjugate_convention="maths") 

1346 1 - 2*I 

1347 >>> q.dot(q1, hermitian=True, conjugate_convention="physics") 

1348 1 + 2*I 

1349 

1350 

1351 See Also 

1352 ======== 

1353 

1354 cross 

1355 multiply 

1356 multiply_elementwise 

1357 """ 

1358 from .dense import Matrix 

1359 

1360 if not isinstance(b, MatrixBase): 

1361 if is_sequence(b): 

1362 if len(b) != self.cols and len(b) != self.rows: 

1363 raise ShapeError( 

1364 "Dimensions incorrect for dot product: %s, %s" % ( 

1365 self.shape, len(b))) 

1366 return self.dot(Matrix(b)) 

1367 else: 

1368 raise TypeError( 

1369 "`b` must be an ordered iterable or Matrix, not %s." % 

1370 type(b)) 

1371 

1372 if (1 not in self.shape) or (1 not in b.shape): 

1373 raise ShapeError 

1374 if len(self) != len(b): 

1375 raise ShapeError( 

1376 "Dimensions incorrect for dot product: %s, %s" % (self.shape, b.shape)) 

1377 

1378 mat = self 

1379 n = len(mat) 

1380 if mat.shape != (1, n): 

1381 mat = mat.reshape(1, n) 

1382 if b.shape != (n, 1): 

1383 b = b.reshape(n, 1) 

1384 

1385 # Now ``mat`` is a row vector and ``b`` is a column vector. 

1386 

1387 # If it so happens that only conjugate_convention is passed 

1388 # then automatically set hermitian to True. If only hermitian 

1389 # is true but no conjugate_convention is not passed then 

1390 # automatically set it to ``"maths"`` 

1391 

1392 if conjugate_convention is not None and hermitian is None: 

1393 hermitian = True 

1394 if hermitian and conjugate_convention is None: 

1395 conjugate_convention = "maths" 

1396 

1397 if hermitian == True: 

1398 if conjugate_convention in ("maths", "left", "math"): 

1399 mat = mat.conjugate() 

1400 elif conjugate_convention in ("physics", "right"): 

1401 b = b.conjugate() 

1402 else: 

1403 raise ValueError("Unknown conjugate_convention was entered." 

1404 " conjugate_convention must be one of the" 

1405 " following: math, maths, left, physics or right.") 

1406 return (mat * b)[0] 

1407 

1408 def dual(self): 

1409 """Returns the dual of a matrix. 

1410 

1411 A dual of a matrix is: 

1412 

1413 ``(1/2)*levicivita(i, j, k, l)*M(k, l)`` summed over indices `k` and `l` 

1414 

1415 Since the levicivita method is anti_symmetric for any pairwise 

1416 exchange of indices, the dual of a symmetric matrix is the zero 

1417 matrix. Strictly speaking the dual defined here assumes that the 

1418 'matrix' `M` is a contravariant anti_symmetric second rank tensor, 

1419 so that the dual is a covariant second rank tensor. 

1420 

1421 """ 

1422 from sympy.matrices import zeros 

1423 

1424 M, n = self[:, :], self.rows 

1425 work = zeros(n) 

1426 if self.is_symmetric(): 

1427 return work 

1428 

1429 for i in range(1, n): 

1430 for j in range(1, n): 

1431 acum = 0 

1432 for k in range(1, n): 

1433 acum += LeviCivita(i, j, 0, k) * M[0, k] 

1434 work[i, j] = acum 

1435 work[j, i] = -acum 

1436 

1437 for l in range(1, n): 

1438 acum = 0 

1439 for a in range(1, n): 

1440 for b in range(1, n): 

1441 acum += LeviCivita(0, l, a, b) * M[a, b] 

1442 acum /= 2 

1443 work[0, l] = -acum 

1444 work[l, 0] = acum 

1445 

1446 return work 

1447 

1448 def _eval_matrix_exp_jblock(self): 

1449 """A helper function to compute an exponential of a Jordan block 

1450 matrix 

1451 

1452 Examples 

1453 ======== 

1454 

1455 >>> from sympy import Symbol, Matrix 

1456 >>> l = Symbol('lamda') 

1457 

1458 A trivial example of 1*1 Jordan block: 

1459 

1460 >>> m = Matrix.jordan_block(1, l) 

1461 >>> m._eval_matrix_exp_jblock() 

1462 Matrix([[exp(lamda)]]) 

1463 

1464 An example of 3*3 Jordan block: 

1465 

1466 >>> m = Matrix.jordan_block(3, l) 

1467 >>> m._eval_matrix_exp_jblock() 

1468 Matrix([ 

1469 [exp(lamda), exp(lamda), exp(lamda)/2], 

1470 [ 0, exp(lamda), exp(lamda)], 

1471 [ 0, 0, exp(lamda)]]) 

1472 

1473 References 

1474 ========== 

1475 

1476 .. [1] https://en.wikipedia.org/wiki/Matrix_function#Jordan_decomposition 

1477 """ 

1478 size = self.rows 

1479 l = self[0, 0] 

1480 exp_l = exp(l) 

1481 

1482 bands = {i: exp_l / factorial(i) for i in range(size)} 

1483 

1484 from .sparsetools import banded 

1485 return self.__class__(banded(size, bands)) 

1486 

1487 

1488 def analytic_func(self, f, x): 

1489 """ 

1490 Computes f(A) where A is a Square Matrix 

1491 and f is an analytic function. 

1492 

1493 Examples 

1494 ======== 

1495 

1496 >>> from sympy import Symbol, Matrix, S, log 

1497 

1498 >>> x = Symbol('x') 

1499 >>> m = Matrix([[S(5)/4, S(3)/4], [S(3)/4, S(5)/4]]) 

1500 >>> f = log(x) 

1501 >>> m.analytic_func(f, x) 

1502 Matrix([ 

1503 [ 0, log(2)], 

1504 [log(2), 0]]) 

1505 

1506 Parameters 

1507 ========== 

1508 

1509 f : Expr 

1510 Analytic Function 

1511 x : Symbol 

1512 parameter of f 

1513 

1514 """ 

1515 

1516 f, x = _sympify(f), _sympify(x) 

1517 if not self.is_square: 

1518 raise NonSquareMatrixError 

1519 if not x.is_symbol: 

1520 raise ValueError("{} must be a symbol.".format(x)) 

1521 if x not in f.free_symbols: 

1522 raise ValueError( 

1523 "{} must be a parameter of {}.".format(x, f)) 

1524 if x in self.free_symbols: 

1525 raise ValueError( 

1526 "{} must not be a parameter of {}.".format(x, self)) 

1527 

1528 eigen = self.eigenvals() 

1529 max_mul = max(eigen.values()) 

1530 derivative = {} 

1531 dd = f 

1532 for i in range(max_mul - 1): 

1533 dd = diff(dd, x) 

1534 derivative[i + 1] = dd 

1535 n = self.shape[0] 

1536 r = self.zeros(n) 

1537 f_val = self.zeros(n, 1) 

1538 row = 0 

1539 

1540 for i in eigen: 

1541 mul = eigen[i] 

1542 f_val[row] = f.subs(x, i) 

1543 if f_val[row].is_number and not f_val[row].is_complex: 

1544 raise ValueError( 

1545 "Cannot evaluate the function because the " 

1546 "function {} is not analytic at the given " 

1547 "eigenvalue {}".format(f, f_val[row])) 

1548 val = 1 

1549 for a in range(n): 

1550 r[row, a] = val 

1551 val *= i 

1552 if mul > 1: 

1553 coe = [1 for ii in range(n)] 

1554 deri = 1 

1555 while mul > 1: 

1556 row = row + 1 

1557 mul -= 1 

1558 d_i = derivative[deri].subs(x, i) 

1559 if d_i.is_number and not d_i.is_complex: 

1560 raise ValueError( 

1561 "Cannot evaluate the function because the " 

1562 "derivative {} is not analytic at the given " 

1563 "eigenvalue {}".format(derivative[deri], d_i)) 

1564 f_val[row] = d_i 

1565 for a in range(n): 

1566 if a - deri + 1 <= 0: 

1567 r[row, a] = 0 

1568 coe[a] = 0 

1569 continue 

1570 coe[a] = coe[a]*(a - deri + 1) 

1571 r[row, a] = coe[a]*pow(i, a - deri) 

1572 deri += 1 

1573 row += 1 

1574 c = r.solve(f_val) 

1575 ans = self.zeros(n) 

1576 pre = self.eye(n) 

1577 for i in range(n): 

1578 ans = ans + c[i]*pre 

1579 pre *= self 

1580 return ans 

1581 

1582 

1583 def exp(self): 

1584 """Return the exponential of a square matrix. 

1585 

1586 Examples 

1587 ======== 

1588 

1589 >>> from sympy import Symbol, Matrix 

1590 

1591 >>> t = Symbol('t') 

1592 >>> m = Matrix([[0, 1], [-1, 0]]) * t 

1593 >>> m.exp() 

1594 Matrix([ 

1595 [ exp(I*t)/2 + exp(-I*t)/2, -I*exp(I*t)/2 + I*exp(-I*t)/2], 

1596 [I*exp(I*t)/2 - I*exp(-I*t)/2, exp(I*t)/2 + exp(-I*t)/2]]) 

1597 """ 

1598 if not self.is_square: 

1599 raise NonSquareMatrixError( 

1600 "Exponentiation is valid only for square matrices") 

1601 try: 

1602 P, J = self.jordan_form() 

1603 cells = J.get_diag_blocks() 

1604 except MatrixError: 

1605 raise NotImplementedError( 

1606 "Exponentiation is implemented only for matrices for which the Jordan normal form can be computed") 

1607 

1608 blocks = [cell._eval_matrix_exp_jblock() for cell in cells] 

1609 from sympy.matrices import diag 

1610 eJ = diag(*blocks) 

1611 # n = self.rows 

1612 ret = P.multiply(eJ, dotprodsimp=None).multiply(P.inv(), dotprodsimp=None) 

1613 if all(value.is_real for value in self.values()): 

1614 return type(self)(re(ret)) 

1615 else: 

1616 return type(self)(ret) 

1617 

1618 def _eval_matrix_log_jblock(self): 

1619 """Helper function to compute logarithm of a jordan block. 

1620 

1621 Examples 

1622 ======== 

1623 

1624 >>> from sympy import Symbol, Matrix 

1625 >>> l = Symbol('lamda') 

1626 

1627 A trivial example of 1*1 Jordan block: 

1628 

1629 >>> m = Matrix.jordan_block(1, l) 

1630 >>> m._eval_matrix_log_jblock() 

1631 Matrix([[log(lamda)]]) 

1632 

1633 An example of 3*3 Jordan block: 

1634 

1635 >>> m = Matrix.jordan_block(3, l) 

1636 >>> m._eval_matrix_log_jblock() 

1637 Matrix([ 

1638 [log(lamda), 1/lamda, -1/(2*lamda**2)], 

1639 [ 0, log(lamda), 1/lamda], 

1640 [ 0, 0, log(lamda)]]) 

1641 """ 

1642 size = self.rows 

1643 l = self[0, 0] 

1644 

1645 if l.is_zero: 

1646 raise MatrixError( 

1647 'Could not take logarithm or reciprocal for the given ' 

1648 'eigenvalue {}'.format(l)) 

1649 

1650 bands = {0: log(l)} 

1651 for i in range(1, size): 

1652 bands[i] = -((-l) ** -i) / i 

1653 

1654 from .sparsetools import banded 

1655 return self.__class__(banded(size, bands)) 

1656 

1657 def log(self, simplify=cancel): 

1658 """Return the logarithm of a square matrix. 

1659 

1660 Parameters 

1661 ========== 

1662 

1663 simplify : function, bool 

1664 The function to simplify the result with. 

1665 

1666 Default is ``cancel``, which is effective to reduce the 

1667 expression growing for taking reciprocals and inverses for 

1668 symbolic matrices. 

1669 

1670 Examples 

1671 ======== 

1672 

1673 >>> from sympy import S, Matrix 

1674 

1675 Examples for positive-definite matrices: 

1676 

1677 >>> m = Matrix([[1, 1], [0, 1]]) 

1678 >>> m.log() 

1679 Matrix([ 

1680 [0, 1], 

1681 [0, 0]]) 

1682 

1683 >>> m = Matrix([[S(5)/4, S(3)/4], [S(3)/4, S(5)/4]]) 

1684 >>> m.log() 

1685 Matrix([ 

1686 [ 0, log(2)], 

1687 [log(2), 0]]) 

1688 

1689 Examples for non positive-definite matrices: 

1690 

1691 >>> m = Matrix([[S(3)/4, S(5)/4], [S(5)/4, S(3)/4]]) 

1692 >>> m.log() 

1693 Matrix([ 

1694 [ I*pi/2, log(2) - I*pi/2], 

1695 [log(2) - I*pi/2, I*pi/2]]) 

1696 

1697 >>> m = Matrix( 

1698 ... [[0, 0, 0, 1], 

1699 ... [0, 0, 1, 0], 

1700 ... [0, 1, 0, 0], 

1701 ... [1, 0, 0, 0]]) 

1702 >>> m.log() 

1703 Matrix([ 

1704 [ I*pi/2, 0, 0, -I*pi/2], 

1705 [ 0, I*pi/2, -I*pi/2, 0], 

1706 [ 0, -I*pi/2, I*pi/2, 0], 

1707 [-I*pi/2, 0, 0, I*pi/2]]) 

1708 """ 

1709 if not self.is_square: 

1710 raise NonSquareMatrixError( 

1711 "Logarithm is valid only for square matrices") 

1712 

1713 try: 

1714 if simplify: 

1715 P, J = simplify(self).jordan_form() 

1716 else: 

1717 P, J = self.jordan_form() 

1718 

1719 cells = J.get_diag_blocks() 

1720 except MatrixError: 

1721 raise NotImplementedError( 

1722 "Logarithm is implemented only for matrices for which " 

1723 "the Jordan normal form can be computed") 

1724 

1725 blocks = [ 

1726 cell._eval_matrix_log_jblock() 

1727 for cell in cells] 

1728 from sympy.matrices import diag 

1729 eJ = diag(*blocks) 

1730 

1731 if simplify: 

1732 ret = simplify(P * eJ * simplify(P.inv())) 

1733 ret = self.__class__(ret) 

1734 else: 

1735 ret = P * eJ * P.inv() 

1736 

1737 return ret 

1738 

1739 def is_nilpotent(self): 

1740 """Checks if a matrix is nilpotent. 

1741 

1742 A matrix B is nilpotent if for some integer k, B**k is 

1743 a zero matrix. 

1744 

1745 Examples 

1746 ======== 

1747 

1748 >>> from sympy import Matrix 

1749 >>> a = Matrix([[0, 0, 0], [1, 0, 0], [1, 1, 0]]) 

1750 >>> a.is_nilpotent() 

1751 True 

1752 

1753 >>> a = Matrix([[1, 0, 1], [1, 0, 0], [1, 1, 0]]) 

1754 >>> a.is_nilpotent() 

1755 False 

1756 """ 

1757 if not self: 

1758 return True 

1759 if not self.is_square: 

1760 raise NonSquareMatrixError( 

1761 "Nilpotency is valid only for square matrices") 

1762 x = uniquely_named_symbol('x', self, modify=lambda s: '_' + s) 

1763 p = self.charpoly(x) 

1764 if p.args[0] == x ** self.rows: 

1765 return True 

1766 return False 

1767 

1768 def key2bounds(self, keys): 

1769 """Converts a key with potentially mixed types of keys (integer and slice) 

1770 into a tuple of ranges and raises an error if any index is out of ``self``'s 

1771 range. 

1772 

1773 See Also 

1774 ======== 

1775 

1776 key2ij 

1777 """ 

1778 islice, jslice = [isinstance(k, slice) for k in keys] 

1779 if islice: 

1780 if not self.rows: 

1781 rlo = rhi = 0 

1782 else: 

1783 rlo, rhi = keys[0].indices(self.rows)[:2] 

1784 else: 

1785 rlo = a2idx(keys[0], self.rows) 

1786 rhi = rlo + 1 

1787 if jslice: 

1788 if not self.cols: 

1789 clo = chi = 0 

1790 else: 

1791 clo, chi = keys[1].indices(self.cols)[:2] 

1792 else: 

1793 clo = a2idx(keys[1], self.cols) 

1794 chi = clo + 1 

1795 return rlo, rhi, clo, chi 

1796 

1797 def key2ij(self, key): 

1798 """Converts key into canonical form, converting integers or indexable 

1799 items into valid integers for ``self``'s range or returning slices 

1800 unchanged. 

1801 

1802 See Also 

1803 ======== 

1804 

1805 key2bounds 

1806 """ 

1807 if is_sequence(key): 

1808 if not len(key) == 2: 

1809 raise TypeError('key must be a sequence of length 2') 

1810 return [a2idx(i, n) if not isinstance(i, slice) else i 

1811 for i, n in zip(key, self.shape)] 

1812 elif isinstance(key, slice): 

1813 return key.indices(len(self))[:2] 

1814 else: 

1815 return divmod(a2idx(key, len(self)), self.cols) 

1816 

1817 def normalized(self, iszerofunc=_iszero): 

1818 """Return the normalized version of ``self``. 

1819 

1820 Parameters 

1821 ========== 

1822 

1823 iszerofunc : Function, optional 

1824 A function to determine whether ``self`` is a zero vector. 

1825 The default ``_iszero`` tests to see if each element is 

1826 exactly zero. 

1827 

1828 Returns 

1829 ======= 

1830 

1831 Matrix 

1832 Normalized vector form of ``self``. 

1833 It has the same length as a unit vector. However, a zero vector 

1834 will be returned for a vector with norm 0. 

1835 

1836 Raises 

1837 ====== 

1838 

1839 ShapeError 

1840 If the matrix is not in a vector form. 

1841 

1842 See Also 

1843 ======== 

1844 

1845 norm 

1846 """ 

1847 if self.rows != 1 and self.cols != 1: 

1848 raise ShapeError("A Matrix must be a vector to normalize.") 

1849 norm = self.norm() 

1850 if iszerofunc(norm): 

1851 out = self.zeros(self.rows, self.cols) 

1852 else: 

1853 out = self.applyfunc(lambda i: i / norm) 

1854 return out 

1855 

1856 def norm(self, ord=None): 

1857 """Return the Norm of a Matrix or Vector. 

1858 

1859 In the simplest case this is the geometric size of the vector 

1860 Other norms can be specified by the ord parameter 

1861 

1862 

1863 ===== ============================ ========================== 

1864 ord norm for matrices norm for vectors 

1865 ===== ============================ ========================== 

1866 None Frobenius norm 2-norm 

1867 'fro' Frobenius norm - does not exist 

1868 inf maximum row sum max(abs(x)) 

1869 -inf -- min(abs(x)) 

1870 1 maximum column sum as below 

1871 -1 -- as below 

1872 2 2-norm (largest sing. value) as below 

1873 -2 smallest singular value as below 

1874 other - does not exist sum(abs(x)**ord)**(1./ord) 

1875 ===== ============================ ========================== 

1876 

1877 Examples 

1878 ======== 

1879 

1880 >>> from sympy import Matrix, Symbol, trigsimp, cos, sin, oo 

1881 >>> x = Symbol('x', real=True) 

1882 >>> v = Matrix([cos(x), sin(x)]) 

1883 >>> trigsimp( v.norm() ) 

1884 1 

1885 >>> v.norm(10) 

1886 (sin(x)**10 + cos(x)**10)**(1/10) 

1887 >>> A = Matrix([[1, 1], [1, 1]]) 

1888 >>> A.norm(1) # maximum sum of absolute values of A is 2 

1889 2 

1890 >>> A.norm(2) # Spectral norm (max of |Ax|/|x| under 2-vector-norm) 

1891 2 

1892 >>> A.norm(-2) # Inverse spectral norm (smallest singular value) 

1893 0 

1894 >>> A.norm() # Frobenius Norm 

1895 2 

1896 >>> A.norm(oo) # Infinity Norm 

1897 2 

1898 >>> Matrix([1, -2]).norm(oo) 

1899 2 

1900 >>> Matrix([-1, 2]).norm(-oo) 

1901 1 

1902 

1903 See Also 

1904 ======== 

1905 

1906 normalized 

1907 """ 

1908 # Row or Column Vector Norms 

1909 vals = list(self.values()) or [0] 

1910 if S.One in self.shape: 

1911 if ord in (2, None): # Common case sqrt(<x, x>) 

1912 return sqrt(Add(*(abs(i) ** 2 for i in vals))) 

1913 

1914 elif ord == 1: # sum(abs(x)) 

1915 return Add(*(abs(i) for i in vals)) 

1916 

1917 elif ord is S.Infinity: # max(abs(x)) 

1918 return Max(*[abs(i) for i in vals]) 

1919 

1920 elif ord is S.NegativeInfinity: # min(abs(x)) 

1921 return Min(*[abs(i) for i in vals]) 

1922 

1923 # Otherwise generalize the 2-norm, Sum(x_i**ord)**(1/ord) 

1924 # Note that while useful this is not mathematically a norm 

1925 try: 

1926 return Pow(Add(*(abs(i) ** ord for i in vals)), S.One / ord) 

1927 except (NotImplementedError, TypeError): 

1928 raise ValueError("Expected order to be Number, Symbol, oo") 

1929 

1930 # Matrix Norms 

1931 else: 

1932 if ord == 1: # Maximum column sum 

1933 m = self.applyfunc(abs) 

1934 return Max(*[sum(m.col(i)) for i in range(m.cols)]) 

1935 

1936 elif ord == 2: # Spectral Norm 

1937 # Maximum singular value 

1938 return Max(*self.singular_values()) 

1939 

1940 elif ord == -2: 

1941 # Minimum singular value 

1942 return Min(*self.singular_values()) 

1943 

1944 elif ord is S.Infinity: # Infinity Norm - Maximum row sum 

1945 m = self.applyfunc(abs) 

1946 return Max(*[sum(m.row(i)) for i in range(m.rows)]) 

1947 

1948 elif (ord is None or isinstance(ord, 

1949 str) and ord.lower() in 

1950 ['f', 'fro', 'frobenius', 'vector']): 

1951 # Reshape as vector and send back to norm function 

1952 return self.vec().norm(ord=2) 

1953 

1954 else: 

1955 raise NotImplementedError("Matrix Norms under development") 

1956 

1957 def print_nonzero(self, symb="X"): 

1958 """Shows location of non-zero entries for fast shape lookup. 

1959 

1960 Examples 

1961 ======== 

1962 

1963 >>> from sympy import Matrix, eye 

1964 >>> m = Matrix(2, 3, lambda i, j: i*3+j) 

1965 >>> m 

1966 Matrix([ 

1967 [0, 1, 2], 

1968 [3, 4, 5]]) 

1969 >>> m.print_nonzero() 

1970 [ XX] 

1971 [XXX] 

1972 >>> m = eye(4) 

1973 >>> m.print_nonzero("x") 

1974 [x ] 

1975 [ x ] 

1976 [ x ] 

1977 [ x] 

1978 

1979 """ 

1980 s = [] 

1981 for i in range(self.rows): 

1982 line = [] 

1983 for j in range(self.cols): 

1984 if self[i, j] == 0: 

1985 line.append(" ") 

1986 else: 

1987 line.append(str(symb)) 

1988 s.append("[%s]" % ''.join(line)) 

1989 print('\n'.join(s)) 

1990 

1991 def project(self, v): 

1992 """Return the projection of ``self`` onto the line containing ``v``. 

1993 

1994 Examples 

1995 ======== 

1996 

1997 >>> from sympy import Matrix, S, sqrt 

1998 >>> V = Matrix([sqrt(3)/2, S.Half]) 

1999 >>> x = Matrix([[1, 0]]) 

2000 >>> V.project(x) 

2001 Matrix([[sqrt(3)/2, 0]]) 

2002 >>> V.project(-x) 

2003 Matrix([[sqrt(3)/2, 0]]) 

2004 """ 

2005 return v * (self.dot(v) / v.dot(v)) 

2006 

2007 def table(self, printer, rowstart='[', rowend=']', rowsep='\n', 

2008 colsep=', ', align='right'): 

2009 r""" 

2010 String form of Matrix as a table. 

2011 

2012 ``printer`` is the printer to use for on the elements (generally 

2013 something like StrPrinter()) 

2014 

2015 ``rowstart`` is the string used to start each row (by default '['). 

2016 

2017 ``rowend`` is the string used to end each row (by default ']'). 

2018 

2019 ``rowsep`` is the string used to separate rows (by default a newline). 

2020 

2021 ``colsep`` is the string used to separate columns (by default ', '). 

2022 

2023 ``align`` defines how the elements are aligned. Must be one of 'left', 

2024 'right', or 'center'. You can also use '<', '>', and '^' to mean the 

2025 same thing, respectively. 

2026 

2027 This is used by the string printer for Matrix. 

2028 

2029 Examples 

2030 ======== 

2031 

2032 >>> from sympy import Matrix, StrPrinter 

2033 >>> M = Matrix([[1, 2], [-33, 4]]) 

2034 >>> printer = StrPrinter() 

2035 >>> M.table(printer) 

2036 '[ 1, 2]\n[-33, 4]' 

2037 >>> print(M.table(printer)) 

2038 [ 1, 2] 

2039 [-33, 4] 

2040 >>> print(M.table(printer, rowsep=',\n')) 

2041 [ 1, 2], 

2042 [-33, 4] 

2043 >>> print('[%s]' % M.table(printer, rowsep=',\n')) 

2044 [[ 1, 2], 

2045 [-33, 4]] 

2046 >>> print(M.table(printer, colsep=' ')) 

2047 [ 1 2] 

2048 [-33 4] 

2049 >>> print(M.table(printer, align='center')) 

2050 [ 1 , 2] 

2051 [-33, 4] 

2052 >>> print(M.table(printer, rowstart='{', rowend='}')) 

2053 { 1, 2} 

2054 {-33, 4} 

2055 """ 

2056 # Handle zero dimensions: 

2057 if S.Zero in self.shape: 

2058 return '[]' 

2059 # Build table of string representations of the elements 

2060 res = [] 

2061 # Track per-column max lengths for pretty alignment 

2062 maxlen = [0] * self.cols 

2063 for i in range(self.rows): 

2064 res.append([]) 

2065 for j in range(self.cols): 

2066 s = printer._print(self[i, j]) 

2067 res[-1].append(s) 

2068 maxlen[j] = max(len(s), maxlen[j]) 

2069 # Patch strings together 

2070 align = { 

2071 'left': 'ljust', 

2072 'right': 'rjust', 

2073 'center': 'center', 

2074 '<': 'ljust', 

2075 '>': 'rjust', 

2076 '^': 'center', 

2077 }[align] 

2078 for i, row in enumerate(res): 

2079 for j, elem in enumerate(row): 

2080 row[j] = getattr(elem, align)(maxlen[j]) 

2081 res[i] = rowstart + colsep.join(row) + rowend 

2082 return rowsep.join(res) 

2083 

2084 def rank_decomposition(self, iszerofunc=_iszero, simplify=False): 

2085 return _rank_decomposition(self, iszerofunc=iszerofunc, 

2086 simplify=simplify) 

2087 

2088 def cholesky(self, hermitian=True): 

2089 raise NotImplementedError('This function is implemented in DenseMatrix or SparseMatrix') 

2090 

2091 def LDLdecomposition(self, hermitian=True): 

2092 raise NotImplementedError('This function is implemented in DenseMatrix or SparseMatrix') 

2093 

2094 def LUdecomposition(self, iszerofunc=_iszero, simpfunc=None, 

2095 rankcheck=False): 

2096 return _LUdecomposition(self, iszerofunc=iszerofunc, simpfunc=simpfunc, 

2097 rankcheck=rankcheck) 

2098 

2099 def LUdecomposition_Simple(self, iszerofunc=_iszero, simpfunc=None, 

2100 rankcheck=False): 

2101 return _LUdecomposition_Simple(self, iszerofunc=iszerofunc, 

2102 simpfunc=simpfunc, rankcheck=rankcheck) 

2103 

2104 def LUdecompositionFF(self): 

2105 return _LUdecompositionFF(self) 

2106 

2107 def singular_value_decomposition(self): 

2108 return _singular_value_decomposition(self) 

2109 

2110 def QRdecomposition(self): 

2111 return _QRdecomposition(self) 

2112 

2113 def upper_hessenberg_decomposition(self): 

2114 return _upper_hessenberg_decomposition(self) 

2115 

2116 def diagonal_solve(self, rhs): 

2117 return _diagonal_solve(self, rhs) 

2118 

2119 def lower_triangular_solve(self, rhs): 

2120 raise NotImplementedError('This function is implemented in DenseMatrix or SparseMatrix') 

2121 

2122 def upper_triangular_solve(self, rhs): 

2123 raise NotImplementedError('This function is implemented in DenseMatrix or SparseMatrix') 

2124 

2125 def cholesky_solve(self, rhs): 

2126 return _cholesky_solve(self, rhs) 

2127 

2128 def LDLsolve(self, rhs): 

2129 return _LDLsolve(self, rhs) 

2130 

2131 def LUsolve(self, rhs, iszerofunc=_iszero): 

2132 return _LUsolve(self, rhs, iszerofunc=iszerofunc) 

2133 

2134 def QRsolve(self, b): 

2135 return _QRsolve(self, b) 

2136 

2137 def gauss_jordan_solve(self, B, freevar=False): 

2138 return _gauss_jordan_solve(self, B, freevar=freevar) 

2139 

2140 def pinv_solve(self, B, arbitrary_matrix=None): 

2141 return _pinv_solve(self, B, arbitrary_matrix=arbitrary_matrix) 

2142 

2143 def solve(self, rhs, method='GJ'): 

2144 return _solve(self, rhs, method=method) 

2145 

2146 def solve_least_squares(self, rhs, method='CH'): 

2147 return _solve_least_squares(self, rhs, method=method) 

2148 

2149 def pinv(self, method='RD'): 

2150 return _pinv(self, method=method) 

2151 

2152 def inv_mod(self, m): 

2153 return _inv_mod(self, m) 

2154 

2155 def inverse_ADJ(self, iszerofunc=_iszero): 

2156 return _inv_ADJ(self, iszerofunc=iszerofunc) 

2157 

2158 def inverse_BLOCK(self, iszerofunc=_iszero): 

2159 return _inv_block(self, iszerofunc=iszerofunc) 

2160 

2161 def inverse_GE(self, iszerofunc=_iszero): 

2162 return _inv_GE(self, iszerofunc=iszerofunc) 

2163 

2164 def inverse_LU(self, iszerofunc=_iszero): 

2165 return _inv_LU(self, iszerofunc=iszerofunc) 

2166 

2167 def inverse_CH(self, iszerofunc=_iszero): 

2168 return _inv_CH(self, iszerofunc=iszerofunc) 

2169 

2170 def inverse_LDL(self, iszerofunc=_iszero): 

2171 return _inv_LDL(self, iszerofunc=iszerofunc) 

2172 

2173 def inverse_QR(self, iszerofunc=_iszero): 

2174 return _inv_QR(self, iszerofunc=iszerofunc) 

2175 

2176 def inv(self, method=None, iszerofunc=_iszero, try_block_diag=False): 

2177 return _inv(self, method=method, iszerofunc=iszerofunc, 

2178 try_block_diag=try_block_diag) 

2179 

2180 def connected_components(self): 

2181 return _connected_components(self) 

2182 

2183 def connected_components_decomposition(self): 

2184 return _connected_components_decomposition(self) 

2185 

2186 def strongly_connected_components(self): 

2187 return _strongly_connected_components(self) 

2188 

2189 def strongly_connected_components_decomposition(self, lower=True): 

2190 return _strongly_connected_components_decomposition(self, lower=lower) 

2191 

2192 _sage_ = Basic._sage_ 

2193 

2194 rank_decomposition.__doc__ = _rank_decomposition.__doc__ 

2195 cholesky.__doc__ = _cholesky.__doc__ 

2196 LDLdecomposition.__doc__ = _LDLdecomposition.__doc__ 

2197 LUdecomposition.__doc__ = _LUdecomposition.__doc__ 

2198 LUdecomposition_Simple.__doc__ = _LUdecomposition_Simple.__doc__ 

2199 LUdecompositionFF.__doc__ = _LUdecompositionFF.__doc__ 

2200 singular_value_decomposition.__doc__ = _singular_value_decomposition.__doc__ 

2201 QRdecomposition.__doc__ = _QRdecomposition.__doc__ 

2202 upper_hessenberg_decomposition.__doc__ = _upper_hessenberg_decomposition.__doc__ 

2203 

2204 diagonal_solve.__doc__ = _diagonal_solve.__doc__ 

2205 lower_triangular_solve.__doc__ = _lower_triangular_solve.__doc__ 

2206 upper_triangular_solve.__doc__ = _upper_triangular_solve.__doc__ 

2207 cholesky_solve.__doc__ = _cholesky_solve.__doc__ 

2208 LDLsolve.__doc__ = _LDLsolve.__doc__ 

2209 LUsolve.__doc__ = _LUsolve.__doc__ 

2210 QRsolve.__doc__ = _QRsolve.__doc__ 

2211 gauss_jordan_solve.__doc__ = _gauss_jordan_solve.__doc__ 

2212 pinv_solve.__doc__ = _pinv_solve.__doc__ 

2213 solve.__doc__ = _solve.__doc__ 

2214 solve_least_squares.__doc__ = _solve_least_squares.__doc__ 

2215 

2216 pinv.__doc__ = _pinv.__doc__ 

2217 inv_mod.__doc__ = _inv_mod.__doc__ 

2218 inverse_ADJ.__doc__ = _inv_ADJ.__doc__ 

2219 inverse_GE.__doc__ = _inv_GE.__doc__ 

2220 inverse_LU.__doc__ = _inv_LU.__doc__ 

2221 inverse_CH.__doc__ = _inv_CH.__doc__ 

2222 inverse_LDL.__doc__ = _inv_LDL.__doc__ 

2223 inverse_QR.__doc__ = _inv_QR.__doc__ 

2224 inverse_BLOCK.__doc__ = _inv_block.__doc__ 

2225 inv.__doc__ = _inv.__doc__ 

2226 

2227 connected_components.__doc__ = _connected_components.__doc__ 

2228 connected_components_decomposition.__doc__ = \ 

2229 _connected_components_decomposition.__doc__ 

2230 strongly_connected_components.__doc__ = \ 

2231 _strongly_connected_components.__doc__ 

2232 strongly_connected_components_decomposition.__doc__ = \ 

2233 _strongly_connected_components_decomposition.__doc__