#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):
Another implementation, based on modifying echelon_form (this isn't PLU,
just LU, and it still has leftover cruft from the echelon_form function
{{{
def _lu_decomposition_in_place_classical(self):
"""
Return a matrix such that the lower-triangular part of the
*pivot columns* is L (diagonal of L is assumed to be all
ones). U is found by taking the upper-triangular (including
diagonal) submatrix of pivot columns, then inserting the
non-pivot columns in the appropriate places. We should have
self=L*U.
EXAMPLES::
sage: t = matrix(QQ, 3, range(9)); t
[0 1 2]
[3 4 5]
[6 7 8]
sage: E = t._echelon_in_place_classical(); t
[ 1 0 -1]
[ 0 1 2]
[ 0 0 0]
"""
tm = verbose('generic in-place Gauss elimination on %s x %s
matrix'%(self._nrows, self._ncols))
cdef Py_ssize_t start_row, c, r, nr, nc, i
if self.fetch('in_echelon_form'):
return
self.check_mutability()
cdef Matrix A
nr = self._nrows
nc = self._ncols
A = self
start_row = 0
pivots = []
for c from 0 <= c < nc:
print "checking column ", c
if PyErr_CheckSignals(): raise KeyboardInterrupt
for r from start_row <= r < nr:
print "checking row ", r
if A.get_unsafe(r, c):
#pivots.append(c)
#a_inverse = ~A.get_unsafe(r,c)
#A.rescale_row(r, a_inverse, c)
# We assume that we do not have to do row swaps
# (later, we'll implement PLU decomposition)
if r != start_row:
raise ValueError, "row reduction required row
swaps, which is not allowed in generic LU decomposition"
#A.swap_rows(r, start_row)
for i from start_row < i < nr:
if A.get_unsafe(i,c):
minus_b = -A.get_unsafe(i,
c)*~A.get_unsafe(r,c)
#print "replacing row ", i, " with i + ",
minus_b, " times row ", start_row, ", starting at column ", c
#print self.str()
A.add_multiple_of_row(i, start_row, minus_b,
c)
A.set_unsafe(i,c,-minus_b)
start_row = start_row + 1
break
#self.cache('pivots', pivots)
#self.cache('in_echelon_form', True)
#self.cache('echelon_form', self)
verbose('done with LU decomposition', tm)
}}}
--
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/3048#comment:4>
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
-~----------~----~----~----~------~----~------~--~---
