> restart;
> with(linalg):
> A := linalg[matrix](3,3, [4,1,1,-1,3,1,-1,1,2]);
> B := vector([9,8,7]);
> for k from 1 by 1 to 3 do
> det(linalg[submatrix](A, 1..k, 1..k));
> end do;
> L1 := linalg[matrix](3,3, [1,0,0,0,1,0,0,0,1]);
> L1[2,1] := -A[2,1]/A[1,1]:
> L1[3,1] := -A[3,1]/A[1,1]:
> L1 = evalm(L1);
> L2 := linalg[matrix](3,3, [1,0,0,0,1,0,0,0,1]);
> A2 := linalg[multiply](L1, A);
> L2[3,2] := -A2[3,2]/A2[2,2]:
> L2 = evalm(L2);
> L = inverse(linalg[multiply](L1, L2));
> 
> A3 := linalg[multiply](L2, A2);
> 
> n := 3;
> LinearAlgebra