我搜索了 Gilbert-Peierls 算法，但没有发现任何有用的东西（好吧，我找到了这个，但它没有按应有的方式工作）。我认为问题出在第二部分，还有那些行：

U(1:k, k) = x(1:k);
L(k:N, k) = x(k:N)/U(k, k);

实际上应该是（根据这个例子）：

U(1:N, k) = x(1:N);
L(k:N, k) = x(k:N)/U(k, k);

另外，在那个例子中，L 是单位矩阵，我觉得有点奇怪。有人可以描述一下算法吗？（有或没有代码）