symfactq(a)

function [f,Q,S] = symfactq(a)
%Symfactq  Symmetric matrix factoring
%
%This program factors a square nxn matrix
%in a product of n+1 symmetric matrices.
%
%usage: [F, Q, S] = symfactq(A)
%
%where F(:,:,1:n-1) will be symmetric orthogonal matrices
%and S will be a symmetric matric "containing" the singular values
%of A such that Q=F(:,:,1)*...*F(:,:,n) and A=Q*S
%
%see also: SVD, EIG, CHOL, PROCRUSTES, private/house

%Paul Godfrey
%Oct. 13, 2006

n=length(a);
E=eye(n);

[u,s,v]=svd(a);
q=u*v';
S=v*s*v';
f(n,n,n)=0;

for k=1:n-1
    x=q(:,1);
    x(1)=x(1)-norm(x);
    p=x'*x;
    if p==0,p=1;end
    h=eye(size(q))-(2/p)*(x*x');
    q=h*q;
    q=q(2:end,2:end);
    H=E;
    H(k:n,k:n)=h;
    f(:,:,k)=H;
end
f(:,:,n)=E;
f(n,n,n)=q;

Q=E;
for k=1:n
    Q=Q*f(:,:,k);
end

return