LDLT pivot

function [ L, d, p ] = LdLT_pivot( A )
  % A is symmetric matrix
  [n, n] = size(A);
  d = [1:n];
  L = eye(n);
  p = [1:n];
  for col = 1:n-1
    display("NEW COL");
    display(A);
    [x, t] = max(abs(diag(A(col:n,col:n))));
    t = t + col - 1;
    display(col);
    display(t);

    p([col t]) = p([t col]);

    A([col t], :) = A([t col], :);
    A(:, [col t]) = A(:, [t col]);

    L([col t], 1:col-1) = L([t col], 1:col-1);

    d(col) = A(col, col);

    L(col+1:n, col) = A(col, col+1:n)' / d(col);
    display("after swap");
    display(A);
    for row = col+1:n
      A(row, row:n) = A(row, row:n) - L(row, col) * A(col, row:n);
    endfor
  endfor
  d(n) = A(n, n);

endfunction