2836.sagews

Sage includes a copy of LAPACK (in the scipy library).  Here we directly call LAPACK's DGETRF Fortran function to compute a PLU decomposition with partial pivoting.  I do this so that you see how this is done with the industry-standard function for computing LU decompositions.

a=matrix(RDF,[[10,2,3],[3,4,6],[2,8,10]])
a

[10.0 2.0 3.0] [ 3.0 4.0 6.0] [ 2.0 8.0 10.0]

import scipy
import scipy.linalg
lu,piv,success=scipy.linalg.flapack.dgetrf(a.rows())

rows=lu.shape[0]
cols=lu.shape[1]
L=copy(matrix(RDF,rows,cols,1)) # need to make copy so I can change it
U=matrix(RDF,rows,cols)
for i in range(rows):
    for j in range(cols):
        if i>j:
            L[i,j]=lu[i,j]
        else:
            U[i,j]=lu[i,j]

P=copy(identity_matrix(rows)) #need to make copy so I can change it
for i,j in enumerate(piv):
    P.swap_rows(i,j)
    
print "P: \n", P
print "\nL: \n",L
print "\nU: \n",U

P: [1 0 0] [0 0 1] [0 1 0] L: [ 1.0 0.0 0.0] [ 0.2 1.0 0.0] [ 0.3 0.447368421053 1.0] U: [ 10.0 2.0 3.0] [ 0.0 7.6 9.4] [ 0.0 0.0 0.894736842105]

P*L*U

[10.0 2.0 3.0] [ 3.0 4.0 6.0] [ 2.0 8.0 10.0]

Sage calls this function (through scipy) behind the scenes when you ask for an LU decomposition of a real double matrix.

p,l,u=a.LU()
p=p.transpose() # Sage uses the book convention, not the convention from LAPACK or our homework assignment

p

[1.0 0.0 0.0] [0.0 0.0 1.0] [0.0 1.0 0.0]

l

[ 1.0 0.0 0.0] [ 0.2 1.0 0.0] [ 0.3 0.447368421053 1.0]

u

[ 10.0 2.0 3.0] [ 0.0 7.6 9.4] [ 0.0 0.0 0.894736842105]