Symmetric pivoting and elimination for a possibly indefinite symmetric matrix using Aasen’s tridiagonalization in Octave - column by column implementation - Problem 21.5 in Trefethen and Bau

AasenByColumn.m

function [L,D,P] = AasenByColumn(A)

	[m,n]=size(A);

	D=tril(A);

	L=eye(m);
	P=eye(m);
	i=1;
	for i=1:m-1

		Li=eye(m);
		Pi=eye(m);

		[sigma, r] = max(abs(D(i+1:m,i)));

		if sigma == 0
			continue
		endif

		r = r + i;
		Pi(:,[i+1,r])=Pi(:,[r,i+1]);

		L=Pi*L*Pi;

		% symmetric row and column swap on the triangular matrix
		if r != i+1
			D([i+1,r],1:i)=D([r,i+1],1:i);
			D(r+1:m,[i+1,r])=D(r+1:m,[r,i+1]);

			tmp=D(r,i+2:r-1);
			D(r,i+2:r-1)=D(i+2:r-1,i+1);
			D(i+2:r-1,i+1)=tmp;

			tmp = D(i+1,i+1);
			D(i+1,i+1) = D(r,r);
			D(r,r) = tmp;
		endif

		L(i+2:m,i+1)=D(i+2:m,i)/D(i+1,i);

		D(i+2:m,i)=0;

		v=zeros(i,1);

		if i > 1
			v(1)=D(1,2)*L(i+1,2);
			for j=2:i-1
				v(j)=D(j,j-1:j)*L(i+1,j-1:j)'+D(j+1,j)*L(i+1,j+1);
			endfor
			v(i)=D(i,i-1:i)*L(i+1,i-1:i)';
		endif

		D(i+1:m,i+1)=D(i+1:m,i+1)-L(i+1:m,1:i)*v-L(i+1:m,i+1)*D(i+1,i)*L(i+1,i)-L(i+1:m,i)*D(i+1,i); 
		
		D(i+2:m,i+1)=D(i+2:m,i+1)-L(i+2:m,i+1)*D(i+1,i+1);

		P=Pi*P;

	endfor

	D=D+tril(D,-1)';

	%max(max(abs(P*A*P'-L*D*L')))
	
endfunction