test_dense.py 9.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350
  1. from sympy.testing.pytest import raises
  2. from sympy.polys import ZZ, QQ
  3. from sympy.polys.matrices.ddm import DDM
  4. from sympy.polys.matrices.dense import (
  5. ddm_transpose,
  6. ddm_iadd, ddm_isub, ddm_ineg, ddm_imatmul, ddm_imul, ddm_irref,
  7. ddm_idet, ddm_iinv, ddm_ilu, ddm_ilu_split, ddm_ilu_solve, ddm_berk)
  8. from sympy.polys.matrices.exceptions import (
  9. DMDomainError,
  10. DMNonInvertibleMatrixError,
  11. DMNonSquareMatrixError,
  12. DMShapeError,
  13. )
  14. def test_ddm_transpose():
  15. a = [[1, 2], [3, 4]]
  16. assert ddm_transpose(a) == [[1, 3], [2, 4]]
  17. def test_ddm_iadd():
  18. a = [[1, 2], [3, 4]]
  19. b = [[5, 6], [7, 8]]
  20. ddm_iadd(a, b)
  21. assert a == [[6, 8], [10, 12]]
  22. def test_ddm_isub():
  23. a = [[1, 2], [3, 4]]
  24. b = [[5, 6], [7, 8]]
  25. ddm_isub(a, b)
  26. assert a == [[-4, -4], [-4, -4]]
  27. def test_ddm_ineg():
  28. a = [[1, 2], [3, 4]]
  29. ddm_ineg(a)
  30. assert a == [[-1, -2], [-3, -4]]
  31. def test_ddm_matmul():
  32. a = [[1, 2], [3, 4]]
  33. ddm_imul(a, 2)
  34. assert a == [[2, 4], [6, 8]]
  35. a = [[1, 2], [3, 4]]
  36. ddm_imul(a, 0)
  37. assert a == [[0, 0], [0, 0]]
  38. def test_ddm_imatmul():
  39. a = [[1, 2, 3], [4, 5, 6]]
  40. b = [[1, 2], [3, 4], [5, 6]]
  41. c1 = [[0, 0], [0, 0]]
  42. ddm_imatmul(c1, a, b)
  43. assert c1 == [[22, 28], [49, 64]]
  44. c2 = [[0, 0, 0], [0, 0, 0], [0, 0, 0]]
  45. ddm_imatmul(c2, b, a)
  46. assert c2 == [[9, 12, 15], [19, 26, 33], [29, 40, 51]]
  47. b3 = [[1], [2], [3]]
  48. c3 = [[0], [0]]
  49. ddm_imatmul(c3, a, b3)
  50. assert c3 == [[14], [32]]
  51. def test_ddm_irref():
  52. # Empty matrix
  53. A = []
  54. Ar = []
  55. pivots = []
  56. assert ddm_irref(A) == pivots
  57. assert A == Ar
  58. # Standard square case
  59. A = [[QQ(0), QQ(1)], [QQ(1), QQ(1)]]
  60. Ar = [[QQ(1), QQ(0)], [QQ(0), QQ(1)]]
  61. pivots = [0, 1]
  62. assert ddm_irref(A) == pivots
  63. assert A == Ar
  64. # m < n case
  65. A = [[QQ(1), QQ(2), QQ(1)], [QQ(3), QQ(4), QQ(1)]]
  66. Ar = [[QQ(1), QQ(0), QQ(-1)], [QQ(0), QQ(1), QQ(1)]]
  67. pivots = [0, 1]
  68. assert ddm_irref(A) == pivots
  69. assert A == Ar
  70. # same m < n but reversed
  71. A = [[QQ(3), QQ(4), QQ(1)], [QQ(1), QQ(2), QQ(1)]]
  72. Ar = [[QQ(1), QQ(0), QQ(-1)], [QQ(0), QQ(1), QQ(1)]]
  73. pivots = [0, 1]
  74. assert ddm_irref(A) == pivots
  75. assert A == Ar
  76. # m > n case
  77. A = [[QQ(1), QQ(0)], [QQ(1), QQ(3)], [QQ(0), QQ(1)]]
  78. Ar = [[QQ(1), QQ(0)], [QQ(0), QQ(1)], [QQ(0), QQ(0)]]
  79. pivots = [0, 1]
  80. assert ddm_irref(A) == pivots
  81. assert A == Ar
  82. # Example with missing pivot
  83. A = [[QQ(1), QQ(0), QQ(1)], [QQ(3), QQ(0), QQ(1)]]
  84. Ar = [[QQ(1), QQ(0), QQ(0)], [QQ(0), QQ(0), QQ(1)]]
  85. pivots = [0, 2]
  86. assert ddm_irref(A) == pivots
  87. assert A == Ar
  88. # Example with missing pivot and no replacement
  89. A = [[QQ(0), QQ(1)], [QQ(0), QQ(2)], [QQ(1), QQ(0)]]
  90. Ar = [[QQ(1), QQ(0)], [QQ(0), QQ(1)], [QQ(0), QQ(0)]]
  91. pivots = [0, 1]
  92. assert ddm_irref(A) == pivots
  93. assert A == Ar
  94. def test_ddm_idet():
  95. A = []
  96. assert ddm_idet(A, ZZ) == ZZ(1)
  97. A = [[ZZ(2)]]
  98. assert ddm_idet(A, ZZ) == ZZ(2)
  99. A = [[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]]
  100. assert ddm_idet(A, ZZ) == ZZ(-2)
  101. A = [[ZZ(1), ZZ(2), ZZ(3)], [ZZ(1), ZZ(2), ZZ(4)], [ZZ(1), ZZ(3), ZZ(5)]]
  102. assert ddm_idet(A, ZZ) == ZZ(-1)
  103. A = [[ZZ(1), ZZ(2), ZZ(3)], [ZZ(1), ZZ(2), ZZ(4)], [ZZ(1), ZZ(2), ZZ(5)]]
  104. assert ddm_idet(A, ZZ) == ZZ(0)
  105. A = [[QQ(1, 2), QQ(1, 2)], [QQ(1, 3), QQ(1, 4)]]
  106. assert ddm_idet(A, QQ) == QQ(-1, 24)
  107. def test_ddm_inv():
  108. A = []
  109. Ainv = []
  110. ddm_iinv(Ainv, A, QQ)
  111. assert Ainv == A
  112. A = []
  113. Ainv = []
  114. raises(DMDomainError, lambda: ddm_iinv(Ainv, A, ZZ))
  115. A = [[QQ(1), QQ(2)]]
  116. Ainv = [[QQ(0), QQ(0)]]
  117. raises(DMNonSquareMatrixError, lambda: ddm_iinv(Ainv, A, QQ))
  118. A = [[QQ(1, 1), QQ(2, 1)], [QQ(3, 1), QQ(4, 1)]]
  119. Ainv = [[QQ(0), QQ(0)], [QQ(0), QQ(0)]]
  120. Ainv_expected = [[QQ(-2, 1), QQ(1, 1)], [QQ(3, 2), QQ(-1, 2)]]
  121. ddm_iinv(Ainv, A, QQ)
  122. assert Ainv == Ainv_expected
  123. A = [[QQ(1, 1), QQ(2, 1)], [QQ(2, 1), QQ(4, 1)]]
  124. Ainv = [[QQ(0), QQ(0)], [QQ(0), QQ(0)]]
  125. raises(DMNonInvertibleMatrixError, lambda: ddm_iinv(Ainv, A, QQ))
  126. def test_ddm_ilu():
  127. A = []
  128. Alu = []
  129. swaps = ddm_ilu(A)
  130. assert A == Alu
  131. assert swaps == []
  132. A = [[]]
  133. Alu = [[]]
  134. swaps = ddm_ilu(A)
  135. assert A == Alu
  136. assert swaps == []
  137. A = [[QQ(1), QQ(2)], [QQ(3), QQ(4)]]
  138. Alu = [[QQ(1), QQ(2)], [QQ(3), QQ(-2)]]
  139. swaps = ddm_ilu(A)
  140. assert A == Alu
  141. assert swaps == []
  142. A = [[QQ(0), QQ(2)], [QQ(3), QQ(4)]]
  143. Alu = [[QQ(3), QQ(4)], [QQ(0), QQ(2)]]
  144. swaps = ddm_ilu(A)
  145. assert A == Alu
  146. assert swaps == [(0, 1)]
  147. A = [[QQ(1), QQ(2), QQ(3)], [QQ(4), QQ(5), QQ(6)], [QQ(7), QQ(8), QQ(9)]]
  148. Alu = [[QQ(1), QQ(2), QQ(3)], [QQ(4), QQ(-3), QQ(-6)], [QQ(7), QQ(2), QQ(0)]]
  149. swaps = ddm_ilu(A)
  150. assert A == Alu
  151. assert swaps == []
  152. A = [[QQ(0), QQ(1), QQ(2)], [QQ(0), QQ(1), QQ(3)], [QQ(1), QQ(1), QQ(2)]]
  153. Alu = [[QQ(1), QQ(1), QQ(2)], [QQ(0), QQ(1), QQ(3)], [QQ(0), QQ(1), QQ(-1)]]
  154. swaps = ddm_ilu(A)
  155. assert A == Alu
  156. assert swaps == [(0, 2)]
  157. A = [[QQ(1), QQ(2), QQ(3)], [QQ(4), QQ(5), QQ(6)]]
  158. Alu = [[QQ(1), QQ(2), QQ(3)], [QQ(4), QQ(-3), QQ(-6)]]
  159. swaps = ddm_ilu(A)
  160. assert A == Alu
  161. assert swaps == []
  162. A = [[QQ(1), QQ(2)], [QQ(3), QQ(4)], [QQ(5), QQ(6)]]
  163. Alu = [[QQ(1), QQ(2)], [QQ(3), QQ(-2)], [QQ(5), QQ(2)]]
  164. swaps = ddm_ilu(A)
  165. assert A == Alu
  166. assert swaps == []
  167. def test_ddm_ilu_split():
  168. U = []
  169. L = []
  170. Uexp = []
  171. Lexp = []
  172. swaps = ddm_ilu_split(L, U, QQ)
  173. assert U == Uexp
  174. assert L == Lexp
  175. assert swaps == []
  176. U = [[]]
  177. L = [[QQ(1)]]
  178. Uexp = [[]]
  179. Lexp = [[QQ(1)]]
  180. swaps = ddm_ilu_split(L, U, QQ)
  181. assert U == Uexp
  182. assert L == Lexp
  183. assert swaps == []
  184. U = [[QQ(1), QQ(2)], [QQ(3), QQ(4)]]
  185. L = [[QQ(1), QQ(0)], [QQ(0), QQ(1)]]
  186. Uexp = [[QQ(1), QQ(2)], [QQ(0), QQ(-2)]]
  187. Lexp = [[QQ(1), QQ(0)], [QQ(3), QQ(1)]]
  188. swaps = ddm_ilu_split(L, U, QQ)
  189. assert U == Uexp
  190. assert L == Lexp
  191. assert swaps == []
  192. U = [[QQ(1), QQ(2), QQ(3)], [QQ(4), QQ(5), QQ(6)]]
  193. L = [[QQ(1), QQ(0)], [QQ(0), QQ(1)]]
  194. Uexp = [[QQ(1), QQ(2), QQ(3)], [QQ(0), QQ(-3), QQ(-6)]]
  195. Lexp = [[QQ(1), QQ(0)], [QQ(4), QQ(1)]]
  196. swaps = ddm_ilu_split(L, U, QQ)
  197. assert U == Uexp
  198. assert L == Lexp
  199. assert swaps == []
  200. U = [[QQ(1), QQ(2)], [QQ(3), QQ(4)], [QQ(5), QQ(6)]]
  201. L = [[QQ(1), QQ(0), QQ(0)], [QQ(0), QQ(1), QQ(0)], [QQ(0), QQ(0), QQ(1)]]
  202. Uexp = [[QQ(1), QQ(2)], [QQ(0), QQ(-2)], [QQ(0), QQ(0)]]
  203. Lexp = [[QQ(1), QQ(0), QQ(0)], [QQ(3), QQ(1), QQ(0)], [QQ(5), QQ(2), QQ(1)]]
  204. swaps = ddm_ilu_split(L, U, QQ)
  205. assert U == Uexp
  206. assert L == Lexp
  207. assert swaps == []
  208. def test_ddm_ilu_solve():
  209. # Basic example
  210. # A = [[QQ(1), QQ(2)], [QQ(3), QQ(4)]]
  211. U = [[QQ(1), QQ(2)], [QQ(0), QQ(-2)]]
  212. L = [[QQ(1), QQ(0)], [QQ(3), QQ(1)]]
  213. swaps = []
  214. b = DDM([[QQ(1)], [QQ(2)]], (2, 1), QQ)
  215. x = DDM([[QQ(0)], [QQ(0)]], (2, 1), QQ)
  216. xexp = DDM([[QQ(0)], [QQ(1, 2)]], (2, 1), QQ)
  217. ddm_ilu_solve(x, L, U, swaps, b)
  218. assert x == xexp
  219. # Example with swaps
  220. # A = [[QQ(0), QQ(2)], [QQ(3), QQ(4)]]
  221. U = [[QQ(3), QQ(4)], [QQ(0), QQ(2)]]
  222. L = [[QQ(1), QQ(0)], [QQ(0), QQ(1)]]
  223. swaps = [(0, 1)]
  224. b = DDM([[QQ(1)], [QQ(2)]], (2, 1), QQ)
  225. x = DDM([[QQ(0)], [QQ(0)]], (2, 1), QQ)
  226. xexp = DDM([[QQ(0)], [QQ(1, 2)]], (2, 1), QQ)
  227. ddm_ilu_solve(x, L, U, swaps, b)
  228. assert x == xexp
  229. # Overdetermined, consistent
  230. # A = DDM([[QQ(1), QQ(2)], [QQ(3), QQ(4)], [QQ(5), QQ(6)]], (3, 2), QQ)
  231. U = [[QQ(1), QQ(2)], [QQ(0), QQ(-2)], [QQ(0), QQ(0)]]
  232. L = [[QQ(1), QQ(0), QQ(0)], [QQ(3), QQ(1), QQ(0)], [QQ(5), QQ(2), QQ(1)]]
  233. swaps = []
  234. b = DDM([[QQ(1)], [QQ(2)], [QQ(3)]], (3, 1), QQ)
  235. x = DDM([[QQ(0)], [QQ(0)]], (2, 1), QQ)
  236. xexp = DDM([[QQ(0)], [QQ(1, 2)]], (2, 1), QQ)
  237. ddm_ilu_solve(x, L, U, swaps, b)
  238. assert x == xexp
  239. # Overdetermined, inconsistent
  240. b = DDM([[QQ(1)], [QQ(2)], [QQ(4)]], (3, 1), QQ)
  241. raises(DMNonInvertibleMatrixError, lambda: ddm_ilu_solve(x, L, U, swaps, b))
  242. # Square, noninvertible
  243. # A = DDM([[QQ(1), QQ(2)], [QQ(1), QQ(2)]], (2, 2), QQ)
  244. U = [[QQ(1), QQ(2)], [QQ(0), QQ(0)]]
  245. L = [[QQ(1), QQ(0)], [QQ(1), QQ(1)]]
  246. swaps = []
  247. b = DDM([[QQ(1)], [QQ(2)]], (2, 1), QQ)
  248. raises(DMNonInvertibleMatrixError, lambda: ddm_ilu_solve(x, L, U, swaps, b))
  249. # Underdetermined
  250. # A = DDM([[QQ(1), QQ(2)]], (1, 2), QQ)
  251. U = [[QQ(1), QQ(2)]]
  252. L = [[QQ(1)]]
  253. swaps = []
  254. b = DDM([[QQ(3)]], (1, 1), QQ)
  255. raises(NotImplementedError, lambda: ddm_ilu_solve(x, L, U, swaps, b))
  256. # Shape mismatch
  257. b3 = DDM([[QQ(1)], [QQ(2)], [QQ(3)]], (3, 1), QQ)
  258. raises(DMShapeError, lambda: ddm_ilu_solve(x, L, U, swaps, b3))
  259. # Empty shape mismatch
  260. U = [[QQ(1)]]
  261. L = [[QQ(1)]]
  262. swaps = []
  263. x = [[QQ(1)]]
  264. b = []
  265. raises(DMShapeError, lambda: ddm_ilu_solve(x, L, U, swaps, b))
  266. # Empty system
  267. U = []
  268. L = []
  269. swaps = []
  270. b = []
  271. x = []
  272. ddm_ilu_solve(x, L, U, swaps, b)
  273. assert x == []
  274. def test_ddm_charpoly():
  275. A = []
  276. assert ddm_berk(A, ZZ) == [[ZZ(1)]]
  277. A = [[ZZ(1), ZZ(2), ZZ(3)], [ZZ(4), ZZ(5), ZZ(6)], [ZZ(7), ZZ(8), ZZ(9)]]
  278. Avec = [[ZZ(1)], [ZZ(-15)], [ZZ(-18)], [ZZ(0)]]
  279. assert ddm_berk(A, ZZ) == Avec
  280. A = DDM([[ZZ(1), ZZ(2)]], (1, 2), ZZ)
  281. raises(DMShapeError, lambda: ddm_berk(A, ZZ))