kill[all];
n:3;
decimales:3;
prec:10^decimales;

A=A:matrix([10,2,-1],
           [-3,0,2],
           [1,1,-5]);
B=B:matrix([27],[-61.5],[-21.5]);
L : ident(n)$;
U : zeromatrix(n,n)$;
for k:1 thru n do (
    U[k,k] : A[k,k]-sum(L[k,r]*U[r,k],r,1,k-1),
    for j:k+1 thru n do (
        U[k,j] : A[k,j]-sum(L[k,r]*U[r,j],r,1,k-1),
        L[j,k] : (A[j,k]-sum(L[j,r]*U[r,k],r,1,k-1))/U[k,k]
    )
) $
print(A, "=", float(round(L*prec)/prec),
     float(round(U*prec)/prec), "=", L.U) $