# 
# Numerical Mathematics and Computing, Fifth Edition
# Ward Cheney & David Kincaid
# Brooks/Cole Publ. Co.
# (c) 2003
#
# file:  lufactor
#
# Solve linear system
#
with(linalg):
A = matrix([[ 3, -13,  9,   3], 
            [-6,   4,  1, -18],
            [ 6,  -2,  2,   4], 
            [12,  -8,  6,  10]]);
b := vector([-19, -34, 16, 26]);
H := augment(A,b);
L := gausselim(A);
x := backsub(L);
LUdecomp(H);

with(LinearAlgebra):
A := Matrix([ [3,-13,9,3],[-6,4,1,-18],
              [6,-2,2,4],[12,-8,6,10]]);
b := Vector([-19,-34,16,26]);
(p, l, u) := LUDecomposition(A);
x := LinearSolve(A, b);