   function [A,TA,trind] = pgenort(m,n)
% generate m random matrices, symm., n by n. 
%  INPUT: m and n
%  the matrices are stored in the rows of the  m by n^2  matrix   A
%  the upper triangular parts only are stored in m by t(n) matrix TA
%   OUTPUT: A
%  Call: function A = pgenort(m,n)
n2=n^2;
rand('seed',sum(100*clock));
[ii,jj]=find(triu(ones(n)));
trind=(jj-1)*n+ii;           % indices for upper triang part of n by n
B=randn(n);B=.5*(B+B');
B=B/norm(B,'fro');
A=reshape(B,1,n2);
TA=B(trind);
TA=TA';
for i=2:m,
   B=randn(n);B=.5*(B+B');
   for k=1:(i-1),
      Bk=reshape(A(k,:),n,n);
      B=B-trace(B*Bk)*Bk;
   end
   B=B/norm(B,'fro');
   A=[A 
   reshape(B,1,n2)];
   TAt=B(trind);
   TAt=TAt';
   TA=[TA
       TAt];
end
%% verify orthogonality of columns
%for i=1:m,
%  for j=1:m,
%    delta(i,j)=trace(reshape(A(i,:),n,n)*reshape(A(j,:),n,n));
%  end
%end
