Merge pull request #355 from isuruf/lu

wrap LUdecomposition

Author: isuruf

symengine/lib/symengine.pxd

@@ -832,7 +832,6 @@ cdef extern from "<symengine/matrix.h>" namespace "SymEngine":
     DenseMatrix* static_cast_DenseMatrix "static_cast<SymEngine::DenseMatrix*>"(const MatrixBase *a)
     void inverse_FFLU "SymEngine::inverse_fraction_free_LU"(const DenseMatrix &A,
         DenseMatrix &B) nogil except +
-    void pivoted_LU (const DenseMatrix &A, DenseMatrix &L, DenseMatrix &U, vector[int] &P) nogil except +
     void pivoted_LU_solve (const DenseMatrix &A, const DenseMatrix &b, DenseMatrix &x) nogil except +
     void inverse_GJ "SymEngine::inverse_gauss_jordan"(const DenseMatrix &A,
         DenseMatrix &B) nogil except +
@@ -857,6 +856,9 @@ cdef extern from "<symengine/matrix.h>" namespace "SymEngine":
     void dot(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &C) nogil
     void cross(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &C) nogil
 
+cdef extern from "<symengine/matrix.h>":
+    void pivoted_LU (const DenseMatrix &A, DenseMatrix &L, DenseMatrix &U, vector[pair[int, int]] &P) nogil except +
+
 cdef extern from "<symengine/ntheory.h>" namespace "SymEngine":
     int probab_prime_p(const Integer &a, int reps)
     RCP[const Integer] nextprime (const Integer &a) nogil

symengine/lib/symengine_wrapper.pyx

@@ -3678,6 +3678,17 @@ cdef class DenseMatrixBase(MatrixBase):
         deref(self.thisptr).LU(deref(L.thisptr), deref(U.thisptr))
         return L, U
 
+    def LUdecomposition(self):
+        cdef DenseMatrixBase L = self.__class__(self.nrows(), self.ncols())
+        cdef DenseMatrixBase U = self.__class__(self.nrows(), self.ncols())
+        cdef vector[pair[int, int]] perm
+        symengine.pivoted_LU(
+            deref(symengine.static_cast_DenseMatrix(self.thisptr)),
+            deref(symengine.static_cast_DenseMatrix(L.thisptr)),
+            deref(symengine.static_cast_DenseMatrix(U.thisptr)),
+            perm)
+        return L, U, perm
+
     def LDL(self):
         cdef DenseMatrixBase L = self.__class__(self.nrows(), self.ncols())
         cdef DenseMatrixBase D = self.__class__(self.nrows(), self.ncols())

symengine/tests/test_matrices.py

@@ -685,3 +685,15 @@ def test_atoms():
     b = Symbol("b")
     X = DenseMatrix([[a, 2], [b, 4]])
     assert X.atoms(Symbol) == set([a, b])
+
+
+def test_LUdecomp():
+    testmat = DenseMatrix([[0, 2, 5, 3],
+                      [3, 3, 7, 4],
+                      [8, 4, 0, 2],
+                      [-2, 6, 3, 4]])
+    L, U, p = testmat.LUdecomposition()
+    res = L*U
+    for orig, new in p:
+        res.row_swap(orig, new)
+    assert res - testmat == zeros(4)