function [Q,D] = CentroSymSchur(A)
% A is an n-by-n centrosymmetric matrix, Q is orthogonal, and
% D = Q'*A*Q is diagonal. If n = 2m, then Q has the form
%
% Q11 Q12
% Q = E = m-by-m exchange permutation
% EQ11 -E*Q12
% If n = 2m+1 then
%
% Q11 Q12 Q13
% Q = Q21 Q22 0 E = m-by-m exchange permutation
% E*Q11 E*Q12 -E*Q13

[n,n] = size(A);
if rem(n,2)==0
m = n/2;
% Structure of A...

% A11 A12
% A =
% E*A12*E E*A11*E
%
% Idea. If (A11 + A12*E)y = lambda* y then

% A11 A12 y y
% = lambda
% E*A12*E E*A11*E Ey Ey

% and if (A11 - A12*E)y = lambda* y then

% A11 A12 y y
% = lambda
% E*A12*E E*A11*E -Ey -Ey

[Q11,D11] = schur(A(1:m,1:m)+A(1:m,n:-1:m+1));
[Q12,D22] = schur(A(1:m,1:m)-A(1:m,n:-1:m+1));
Q = [Q11 Q12; Q11(m:-1:1,:) -Q12(m:-1:1,:)]/sqrt(2);
D = [D11 zeros(m,m); zeros(m,m) D22];
else
m = floor(n/2);
% Structure of A...

% A11 v A12
% A = v' alfa v'*E
% E*A12*E Ev E*A11*E
%
% Idea. If
% (A11+ A12*E) sqrt(2)v y y
% = lambda
% sqrt(2)v' alfa mu mu
% then
%
% A11 v A12 y y
% v' alfa v'*E sqrt(2)mu = lambda sqrt(2)mu
% E*A12*E Ev E*A11*E Ey Ey

% and if (A11 - A12*E)y = lambda*y then

% A11 v A12 y y
% v' alfa v'*E 0 = lambda 0
% E*A12*E Ev E*A11*E -Ey -Ey

[U,D11] = schur(...
[A(1:m,1:m)+A(1:m,n:-1:m+2) sqrt(2)*A(1:m,m+1);...
sqrt(2)*A(m+1,1:m) A(m+1,m+1) ]...
);
[Q13,D33] = schur(A(1:m,1:m)-A(1:m,n:-1:m+2));
% Construct Q
Q = zeros(n,n);
% left half...
U(m+1,:) = U(m+1,:)*sqrt(2);
Q(1:m+1,1:m+1) = U;
Q(m+2:n,1:m+1) = U(m:-1:1,:);
% right half...
Q(1:m,m+2:n) = Q13;
Q(m+2:n,m+2:n) = -Q13(m:-1:1,:);
% normalize...
for k=1:n
Q(:,k) = Q(:,k)/norm(Q(:,k));
end
D = [D11 zeros(m+1,m); zeros(m,m+1) D33];
end