/* * Calculates the LU-decomposition of a square matrix with pivoting. * * input: A square matrix (over any ring) * output: [P,L,U] list */ Define PivotLU(A) U := A; L := Identity(NumRows(A)); P := L; For K := 1 To NumRows(A) - 1 Do PivotFound := FALSE; For I := K To NumRows(A) Do If U[I][K] <> 0 Then Swap := U[K]; U[K] := U[I]; U[I] := Swap; PivotFound := TRUE; Swap := P[K]; P[K] := P[I]; P[I] := Swap; For J := 1 To K - 1 Do Swap := L[K][J]; L[K][J] := L[I][J]; L[I][J] := Swap; EndFor; Break; EndIf; EndFor; If PivotFound Then For J := K + 1 To NumRows(A) Do L[J][K] := U[J][K]/U[K][K]; For I := K To NumRows(A) Do U[J][I] := U[J][I] - L[J][K]*U[K][I]; EndFor; EndFor; EndIf; EndFor; Return [P,L,U]; EndDefine; /* * Creates the elementary matrix, which swaps row i with row j. * * input: M,I,J positive integers * output: T matrix */ Define SwapMatrix(M,I,J) T := Identity(M); T[I][I] := 0; T[J][J] := 0; T[I][J] := 1; T[J][I] := 1; Return T; EndDefine;