1"""
2Solution of equations using dense matrices.
3
4The dense matrix is stored as a list of lists.
5
6"""
7
8import copy
9
10from sympy.core.power import isqrt
11from sympy.core.symbol import symbols
12from sympy.matrices.densetools import (
13    augment, col, conjugate_transpose, eye, rowadd, rowmul)
14from sympy.utilities.exceptions import SymPyDeprecationWarning
15
16SymPyDeprecationWarning(
17    feature="densesolve",
18    issue=12695,
19    deprecated_since_version="1.1").warn()
20
21
22def row_echelon(matlist, K):
23    """
24    Returns the row echelon form of a matrix with diagonal elements
25    reduced to 1.
26
27    Examples
28    ========
29
30    >>> from sympy.matrices.densesolve import row_echelon
31    >>> from sympy import QQ
32    >>> a = [
33    ... [QQ(3), QQ(7), QQ(4)],
34    ... [QQ(2), QQ(4), QQ(5)],
35    ... [QQ(6), QQ(2), QQ(3)]]
36    >>> row_echelon(a, QQ)
37    [[1, 7/3, 4/3], [0, 1, -7/2], [0, 0, 1]]
38
39    See Also
40    ========
41
42    rref
43    """
44    result_matlist = copy.deepcopy(matlist)
45    nrow = len(result_matlist)
46    for i in range(nrow):
47        if (result_matlist[i][i] != 1 and result_matlist[i][i] != 0):
48            rowmul(result_matlist, i, 1/result_matlist[i][i], K)
49        for j in range(i + 1, nrow):
50            if (result_matlist[j][i] != 0):
51                rowadd(result_matlist, j, i, -result_matlist[j][i], K)
52    return result_matlist
53
54
55def rref(matlist, K):
56    """
57    Returns the reduced row echelon form of a Matrix.
58
59    Examples
60    ========
61
62    >>> from sympy.matrices.densesolve import rref
63    >>> from sympy import QQ
64    >>> a = [
65    ... [QQ(1), QQ(2), QQ(1)],
66    ... [QQ(-2), QQ(-3), QQ(1)],
67    ... [QQ(3), QQ(5), QQ(0)]]
68    >>> rref(a, QQ)
69    [[1, 0, -5], [0, 1, 3], [0, 0, 0]]
70
71    See Also
72    ========
73
74    row_echelon
75    """
76    result_matlist = copy.deepcopy(matlist)
77    result_matlist = row_echelon(result_matlist, K)
78    nrow = len(result_matlist)
79    for i in range(nrow):
80        if result_matlist[i][i] == 1:
81            for j in range(i):
82                rowadd(result_matlist, j, i, -result_matlist[j][i], K)
83    return result_matlist
84
85
86def LU(matlist, K, reverse = 0):
87    """
88    It computes the LU decomposition of a matrix and returns L and U
89    matrices.
90
91    Examples
92    ========
93
94    >>> from sympy.matrices.densesolve import LU
95    >>> from sympy import QQ
96    >>> a = [
97    ... [QQ(1), QQ(2), QQ(3)],
98    ... [QQ(2), QQ(-4), QQ(6)],
99    ... [QQ(3), QQ(-9), QQ(-3)]]
100    >>> LU(a, QQ)
101    ([[1, 0, 0], [2, 1, 0], [3, 15/8, 1]], [[1, 2, 3], [0, -8, 0], [0, 0, -12]])
102
103    See Also
104    ========
105
106    upper_triangle
107    lower_triangle
108    """
109    nrow = len(matlist)
110    new_matlist1, new_matlist2 = eye(nrow, K), copy.deepcopy(matlist)
111    for i in range(nrow):
112        for j in range(i + 1, nrow):
113            if (new_matlist2[j][i] != 0):
114                new_matlist1[j][i] = new_matlist2[j][i]/new_matlist2[i][i]
115                rowadd(new_matlist2, j, i, -new_matlist2[j][i]/new_matlist2[i][i], K)
116    return new_matlist1, new_matlist2
117
118
119def cholesky(matlist, K):
120    """
121    Performs the cholesky decomposition of a Hermitian matrix and
122    returns L and it's conjugate transpose.
123
124    Examples
125    ========
126
127    >>> from sympy.matrices.densesolve import cholesky
128    >>> from sympy import QQ
129    >>> cholesky([[QQ(25), QQ(15), QQ(-5)], [QQ(15), QQ(18), QQ(0)], [QQ(-5), QQ(0), QQ(11)]], QQ)
130    ([[5, 0, 0], [3, 3, 0], [-1, 1, 3]], [[5, 3, -1], [0, 3, 1], [0, 0, 3]])
131
132    See Also
133    ========
134
135    cholesky_solve
136    """
137    new_matlist = copy.deepcopy(matlist)
138    nrow = len(new_matlist)
139    L = eye(nrow, K)
140    for i in range(nrow):
141        for j in range(i + 1):
142            a = K.zero
143            for k in range(j):
144                a += L[i][k]*L[j][k]
145            if i == j:
146                L[i][j] = isqrt(new_matlist[i][j] - a)
147            else:
148                L[i][j] = (new_matlist[i][j] - a)/L[j][j]
149    return L, conjugate_transpose(L, K)
150
151
152def LDL(matlist, K):
153    """
154    Performs the LDL decomposition of a hermitian matrix and returns L, D and
155    transpose of L. Only applicable to rational entries.
156
157    Examples
158    ========
159
160    >>> from sympy.matrices.densesolve import LDL
161    >>> from sympy import QQ
162
163    >>> a = [
164    ... [QQ(4), QQ(12), QQ(-16)],
165    ... [QQ(12), QQ(37), QQ(-43)],
166    ... [QQ(-16), QQ(-43), QQ(98)]]
167    >>> LDL(a, QQ)
168    ([[1, 0, 0], [3, 1, 0], [-4, 5, 1]], [[4, 0, 0], [0, 1, 0], [0, 0, 9]], [[1, 3, -4], [0, 1, 5], [0, 0, 1]])
169
170    """
171    new_matlist = copy.deepcopy(matlist)
172    nrow = len(new_matlist)
173    L, D = eye(nrow, K), eye(nrow, K)
174    for i in range(nrow):
175        for j in range(i + 1):
176            a = K.zero
177            for k in range(j):
178                a += L[i][k]*L[j][k]*D[k][k]
179            if i == j:
180                D[j][j] = new_matlist[j][j] - a
181            else:
182                L[i][j] = (new_matlist[i][j] - a)/D[j][j]
183    return L, D, conjugate_transpose(L, K)
184
185
186def upper_triangle(matlist, K):
187    """
188    Transforms a given matrix to an upper triangle matrix by performing
189    row operations on it.
190
191    Examples
192    ========
193
194    >>> from sympy.matrices.densesolve import upper_triangle
195    >>> from sympy import QQ
196    >>> a = [
197    ... [QQ(4,1), QQ(12,1), QQ(-16,1)],
198    ... [QQ(12,1), QQ(37,1), QQ(-43,1)],
199    ... [QQ(-16,1), QQ(-43,1), QQ(98,1)]]
200    >>> upper_triangle(a, QQ)
201    [[4, 12, -16], [0, 1, 5], [0, 0, 9]]
202
203    See Also
204    ========
205
206    LU
207    """
208    copy_matlist = copy.deepcopy(matlist)
209    lower_triangle, upper_triangle = LU(copy_matlist, K)
210    return upper_triangle
211
212
213def lower_triangle(matlist, K):
214    """
215    Transforms a given matrix to a lower triangle matrix by performing
216    row operations on it.
217
218    Examples
219    ========
220
221    >>> from sympy.matrices.densesolve import lower_triangle
222    >>> from sympy import QQ
223    >>> a = [
224    ... [QQ(4,1), QQ(12,1), QQ(-16)],
225    ... [QQ(12,1), QQ(37,1), QQ(-43,1)],
226    ... [QQ(-16,1), QQ(-43,1), QQ(98,1)]]
227    >>> lower_triangle(a, QQ)
228    [[1, 0, 0], [3, 1, 0], [-4, 5, 1]]
229
230    See Also
231    ========
232
233    LU
234    """
235    copy_matlist = copy.deepcopy(matlist)
236    lower_triangle, upper_triangle = LU(copy_matlist, K, reverse = 1)
237    return lower_triangle
238
239
240def rref_solve(matlist, variable, constant, K):
241    """
242    Solves a system of equations using reduced row echelon form given
243    a matrix of coefficients, a vector of variables and a vector of constants.
244
245    Examples
246    ========
247
248    >>> from sympy.matrices.densesolve import rref_solve
249    >>> from sympy import QQ
250    >>> from sympy import Dummy
251    >>> x, y, z = Dummy('x'), Dummy('y'), Dummy('z')
252    >>> coefficients = [
253    ... [QQ(25), QQ(15), QQ(-5)],
254    ... [QQ(15), QQ(18), QQ(0)],
255    ... [QQ(-5), QQ(0), QQ(11)]]
256    >>> constants = [
257    ... [QQ(2)],
258    ... [QQ(3)],
259    ... [QQ(1)]]
260    >>> variables = [
261    ... [x],
262    ... [y],
263    ... [z]]
264    >>> rref_solve(coefficients, variables, constants, QQ)
265    [[-1/225], [23/135], [4/45]]
266
267    See Also
268    ========
269
270    row_echelon
271    augment
272    """
273    new_matlist = copy.deepcopy(matlist)
274    augmented = augment(new_matlist, constant, K)
275    solution = rref(augmented, K)
276    return col(solution, -1)
277
278
279def LU_solve(matlist, variable, constant, K):
280    """
281    Solves a system of equations using  LU decomposition given a matrix
282    of coefficients, a vector of variables and a vector of constants.
283
284    Examples
285    ========
286
287    >>> from sympy.matrices.densesolve import LU_solve
288    >>> from sympy import QQ
289    >>> from sympy import Dummy
290    >>> x, y, z = Dummy('x'), Dummy('y'), Dummy('z')
291    >>> coefficients = [
292    ... [QQ(2), QQ(-1), QQ(-2)],
293    ... [QQ(-4), QQ(6), QQ(3)],
294    ... [QQ(-4), QQ(-2), QQ(8)]]
295    >>> variables = [
296    ... [x],
297    ... [y],
298    ... [z]]
299    >>> constants = [
300    ... [QQ(-1)],
301    ... [QQ(13)],
302    ... [QQ(-6)]]
303    >>> LU_solve(coefficients, variables, constants, QQ)
304    [[2], [3], [1]]
305
306    See Also
307    ========
308
309    LU
310    forward_substitution
311    backward_substitution
312    """
313    new_matlist = copy.deepcopy(matlist)
314    nrow = len(new_matlist)
315    L, U = LU(new_matlist, K)
316    y = [[i] for i in symbols('y:%i' % nrow)]
317    forward_substitution(L, y, constant, K)
318    backward_substitution(U, variable, y, K)
319    return variable
320
321
322def cholesky_solve(matlist, variable, constant, K):
323    """
324    Solves a system of equations using Cholesky decomposition given
325    a matrix of coefficients, a vector of variables and a vector of constants.
326
327    Examples
328    ========
329
330    >>> from sympy.matrices.densesolve import cholesky_solve
331    >>> from sympy import QQ
332    >>> from sympy import Dummy
333    >>> x, y, z = Dummy('x'), Dummy('y'), Dummy('z')
334    >>> coefficients = [
335    ... [QQ(25), QQ(15), QQ(-5)],
336    ... [QQ(15), QQ(18), QQ(0)],
337    ... [QQ(-5), QQ(0), QQ(11)]]
338    >>> variables = [
339    ... [x],
340    ... [y],
341    ... [z]]
342    >>> coefficients = [
343    ... [QQ(2)],
344    ... [QQ(3)],
345    ... [QQ(1)]]
346    >>> cholesky_solve([[QQ(25), QQ(15), QQ(-5)], [QQ(15), QQ(18), QQ(0)], [QQ(-5), QQ(0), QQ(11)]], [[x], [y], [z]], [[QQ(2)], [QQ(3)], [QQ(1)]], QQ)
347    [[-1/225], [23/135], [4/45]]
348
349    See Also
350    ========
351
352    cholesky
353    forward_substitution
354    backward_substitution
355    """
356    new_matlist = copy.deepcopy(matlist)
357    nrow = len(new_matlist)
358    L, Lstar = cholesky(new_matlist, K)
359    y = [[i] for i in symbols('y:%i' % nrow)]
360    forward_substitution(L, y, constant, K)
361    backward_substitution(Lstar, variable, y, K)
362    return variable
363
364
365def forward_substitution(lower_triangle, variable, constant, K):
366    """
367    Performs forward substitution given a lower triangular matrix, a
368    vector of variables and a vector of constants.
369
370    Examples
371    ========
372
373    >>> from sympy.matrices.densesolve import forward_substitution
374    >>> from sympy import QQ
375    >>> from sympy import Dummy
376    >>> x, y, z = Dummy('x'), Dummy('y'), Dummy('z')
377    >>> a = [
378    ... [QQ(1), QQ(0), QQ(0)],
379    ... [QQ(-2), QQ(1), QQ(0)],
380    ... [QQ(-2), QQ(-1), QQ(1)]]
381    >>> variables = [
382    ... [x],
383    ... [y],
384    ... [z]]
385    >>> constants = [
386    ... [QQ(-1)],
387    ... [QQ(13)],
388    ... [QQ(-6)]]
389    >>> forward_substitution(a, variables, constants, QQ)
390    [[-1], [11], [3]]
391
392    See Also
393    ========
394
395    LU_solve
396    cholesky_solve
397    """
398    copy_lower_triangle = copy.deepcopy(lower_triangle)
399    nrow = len(copy_lower_triangle)
400    for i in range(nrow):
401        a = K.zero
402        for j in range(i):
403            a += copy_lower_triangle[i][j]*variable[j][0]
404        variable[i][0] = (constant[i][0] - a)/copy_lower_triangle[i][i]
405    return variable
406
407
408def backward_substitution(upper_triangle, variable, constant, K):
409    """
410    Performs forward substitution given a lower triangular matrix,
411    a vector of variables and a vector constants.
412
413    Examples
414    ========
415
416    >>> from sympy.matrices.densesolve import backward_substitution
417    >>> from sympy import QQ
418    >>> from sympy import Dummy
419    >>> x, y, z = Dummy('x'), Dummy('y'), Dummy('z')
420    >>> a = [
421    ... [QQ(2), QQ(-1), QQ(-2)],
422    ... [QQ(0), QQ(4), QQ(-1)],
423    ... [QQ(0), QQ(0), QQ(3)]]
424    >>> variables = [
425    ... [x],
426    ... [y],
427    ... [z]]
428    >>> constants = [
429    ... [QQ(-1)],
430    ... [QQ(11)],
431    ... [QQ(3)]]
432    >>> backward_substitution(a, variables, constants, QQ)
433    [[2], [3], [1]]
434
435    See Also
436    ========
437
438    LU_solve
439    cholesky_solve
440    """
441    copy_upper_triangle = copy.deepcopy(upper_triangle)
442    nrow = len(copy_upper_triangle)
443    for i in reversed(range(nrow)):
444        a = K.zero
445        for j in reversed(range(i + 1, nrow)):
446            a += copy_upper_triangle[i][j]*variable[j][0]
447        variable[i][0] = (constant[i][0] - a)/copy_upper_triangle[i][i]
448    return variable
449