#3048: [with patch, reviewer requests further work] add generic LU decomposition
----------------------------+-----------------------------------------------
Reporter: mhansen | Owner: mhansen
Type: enhancement | Status: assigned
Priority: minor | Milestone: sage-4.0.1
Component: linear algebra | Keywords:
----------------------------+-----------------------------------------------
Comment(by jason):
Here is another try. Next step is maybe to compare the speed of Mike's
implementation and my reworking of echelon_form.
Then maybe we can change the generic determinant algorithm to use this and
be really fast.
{{{
def _lu_decomposition_(self):
"""
Return the PLU decomposition P, L, and U such that self=P*L*U.
If self is an mxn matrix, then P is an mxm permutation matrix,
L is a mxm unit lower-triangular matrix, and U is an mxn
upper-triangular matrix.
The decomposition is done over the fraction field of
self.base_ring().
EXAMPLES::
"""
tm = verbose('generic in-place LU decomposition on %s x %s
matrix'%(self._nrows, self._ncols))
cdef Py_ssize_t start_row, c, r, nr, nc, i
x = self.fetch('PLU')
if x is not None:
return x[0].matrix(), x[1], x[2]
cdef Matrix A, L
nr = self._nrows
nc = self._ncols
R = self.base_ring().fraction_field()
A = self.change_ring(R).copy()
L = A.new_matrix(nr, nr)
one = R(1)
for r from 0<= r < nr:
L.set_unsafe(r,r,one)
from sage.groups.all import SymmetricGroup
S = SymmetricGroup(nr)
p = S(1)
start_row = 0
pivots = []
for c from 0 <= c < nc:
if PyErr_CheckSignals(): raise KeyboardInterrupt
for r from start_row <= r < nr:
if A.get_unsafe(r, c):
pivots.append(c)
if r != start_row:
A.swap_rows(r, start_row)
p *= S((r+1, start_row+1))
a_inverse = ~A.get_unsafe(start_row,c)
for i from start_row < i < nr:
if A.get_unsafe(i,c):
b = A.get_unsafe(i, c)*a_inverse
A.add_multiple_of_row(i, start_row, -b, c)
L.set_unsafe(i,start_row,b)
start_row += 1
break
self.cache('pivots', pivots)
self.cache('PLU', (p,L,A))
verbose('done with LU decomposition', tm)
return p.matrix(), L, A
}}}
--
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/3048#comment:5>
Sage <http://sagemath.org/>
Sage: Creating a Viable Open Source Alternative to Magma, Maple, Mathematica,
and MATLAB
--~--~---------~--~----~------------~-------~--~----~
You received this message because you are subscribed to the Google Groups
"sage-trac" group.
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to
[email protected]
For more options, visit this group at
http://groups.google.com/group/sage-trac?hl=en
-~----------~----~----~----~------~----~------~--~---
