︠a9fe3b00-7d70-4389-b375-fa36ac6ccb5bi︠
%html
<p>Sage includes a copy of LAPACK (in the scipy library).&nbsp; Here we directly call LAPACK's <a href="http://www.netlib.org/lapack/double/dgetrf.f">DGETRF</a> Fortran function to compute a PLU decomposition with partial pivoting.&nbsp; I do this so that you see how this is done with the industry-standard function for computing LU decompositions.</p>

︡12b2b57a-0849-4eed-b5dc-8943bc96a470︡{"html": "<p>Sage includes a copy of LAPACK (in the scipy library).&nbsp; Here we directly call LAPACK's <a href=\"http://www.netlib.org/lapack/double/dgetrf.f\">DGETRF</a> Fortran function to compute a PLU decomposition with partial pivoting.&nbsp; I do this so that you see how this is done with the industry-standard function for computing LU decompositions.</p>"}︡
︠42a2ceb9-ea76-4f65-b881-23336dae0705︠
a=matrix(RDF,[[10,2,3],[3,4,6],[2,8,10]])
a
︡9468922f-60e1-4a39-845b-5ababd845da4︡{"stdout": "[10.0  2.0  3.0]\n[ 3.0  4.0  6.0]\n[ 2.0  8.0 10.0]"}︡
︠32e2d508-7fca-4ba7-97fe-4d543d285a8c︠
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
︡5d3773d8-5da5-4eeb-83da-484997633d29︡{"stdout": "P: \n[1 0 0]\n[0 0 1]\n[0 1 0]\n\nL: \n[           1.0            0.0            0.0]\n[           0.2            1.0            0.0]\n[           0.3 0.447368421053            1.0]\n\nU: \n[          10.0            2.0            3.0]\n[           0.0            7.6            9.4]\n[           0.0            0.0 0.894736842105]"}︡
︠45c30d02-9362-4cc8-96c9-c2290fb55d1a︠
P*L*U
︡07546de3-384e-49b1-96e4-f8637343c115︡{"stdout": "[10.0  2.0  3.0]\n[ 3.0  4.0  6.0]\n[ 2.0  8.0 10.0]"}︡
︠e00fe280-fa9a-4cdf-9d9b-7a28131e0130i︠
%html
<p>Sage calls this function (through scipy) behind the scenes when you ask for an LU decomposition of a real double matrix.</p>

︡4adf8cb1-b14e-4698-8faf-bd027c7b46ef︡{"html": "<p>Sage calls this function (through scipy) behind the scenes when you ask for an LU decomposition of a real double matrix.</p>"}︡
︠8645a43e-31b8-418f-916f-a8e4f14e2e4d︠
p,l,u=a.LU()
p=p.transpose() # Sage uses the book convention, not the convention from LAPACK or our homework assignment
︡14209739-4981-4c45-99cc-aca2b1efcc1d︡︡
︠9f3715e9-3ebf-4ada-9437-9ef0eba2614f︠
p
︡f8bea2ea-c85c-4dc0-bb46-0f62d1614eea︡{"stdout": "[1.0 0.0 0.0]\n[0.0 0.0 1.0]\n[0.0 1.0 0.0]"}︡
︠df15e81c-06f4-4778-91b8-983984d4c96d︠
l
︡56ffee1e-619c-467c-8ed8-a497d74870a9︡{"stdout": "[           1.0            0.0            0.0]\n[           0.2            1.0            0.0]\n[           0.3 0.447368421053            1.0]"}︡
︠cd46f548-a4a7-4640-8d87-c0106a76515d︠
u
︡c51819d0-ec4d-4c28-b0e7-4954a8fcb859︡{"stdout": "[          10.0            2.0            3.0]\n[           0.0            7.6            9.4]\n[           0.0            0.0 0.894736842105]"}︡