  function x = StructuredSolve(A,b)
% A is a 2n-by-2n symmetric positive definite matrix that satisfies EAE = A
% where E is the 2n-by-2n exchange matrix. 
% b is a 2n-by-1 vector that satisfies b = Eb.
% x solves Ax = b
  
% Let us start by looking what EAE = A means for a 6-by-6 symmetric A:
%
%                 a  b  c    d  e  f
%                 b  g  h    k  m  e
%                 c  h  n    p  k  d
%           A  =
%                 d  k  p    n  h  c
%                 e  m  k    h  g  b
%                 f  e  d    c  b  a
%
% (i) E*matrix*E reverses the order of the matrix rows and columns
% (ii) E*matrix*E = matrix transposed means the matrix is symmetric about
% its antiddiagonal, i.e., the "line" from the (n,1) to (1,n) entry.
% Thus

%             A11     A12
%        A =                    and A21 = E*A12*E        
%             A21   E*A11*E

% where E is the n-by-n exchange matrix.
%
%           b1
% With b = 
%           E*b1
%
% we discover from more small n experimentation that 
%           x1
%      x = 
%           E*x1

% Writing Ax = b as a 2-0by-2 system and equating entries we conclude that
% (A11 + A12*E)x1 = b1 so...
  m = length(b);
  n = m/2;
  x1 = (A(1:n,1:n) + A(1:n,m:-1:n+1))\b(1:n);    % (A11 + A12*E)x1 = b1
  x = [x1 ; x1(n:-1:1)];
  
% -1 if you do this which multiplies work by a factor of three or so:

% E = eye(n,n); E = E(:,n:-1:1); x1 = (A(1:n,1:n) + A(1:n,n+1:m)*E)\b(1:n)