function [U]=conjorth(W,A)
 %[U]=conjorth(W,A)
 % %INPUT
 % W: A linear independet set of column vectors
 % A: A full rank matrix
 %
 %OUTPUT
 %U: A set of mutually A-conjuate vectors that span the same space as W.
 %   If A if the identity matrix, so that A=eye(size(W)), U is a set of
 %   orthogonal vectors and the method is identical to the Gram-Schmidt
 %   orthogonalization.
 %
 %DESCRIPTION OF A-CONJUGACY
 %Two vectors u and v are A-conjugate if u'*(A*v) = 0.  In a more general
 %sense, two elements of a vector space are A-conjugate if <u,Av=0.
 %A-conjugacy arises out of generalized eigenvalue problems of the form:
 %
 % Lu = lambda*A*u
 %
 % If A is self-adjoint and positive definite, and L is self-adjoint, then
 % the eigenvectors  of operator L are A-conjugate to each other.  Such
 % problems arise all the time in vibrational mechanics.  The Ritz method is
 % a common aplication of A-conjugacy of eigenvectors.  In seismology,
 % non-spherical bodies will in general, have their wave-equation
 % eigenvectors A-conjugate to each other, rather than orthogonal, though
 % the elements of A are likely to be quite small for the Earth.
 %
 %USE WITH:
 %Rayleigh-Ritz algorithms, polarization algorithms, Gram-Schmidt
 %orthogonalization, eigenvalue problems.

 [Mw,Nw]=size(W);
 [Ma,Na]=size(A);

 if(Mw<Ma)
    disp(sprintf('U will not span R^(%s) due to size(W)',num2str(Mw)));
 end;

 % MAKING w1 and w2 A-conjugate to each other: Illustrating the method via
 % construction of the first two conjugate-orthogonal vectors:

 % Select a vector from W to be the first element of the A-conjugate set:
 v=W(:,1)/norm(W(:,1));
 % Normalize A*v for use in the projector matrix P that projects onto the
 % range of u
 u=(A*v)/norm(A*v);
 P=(u*u');
 I=eye(size(A));
 % Construct the complementary projector I-P which projects onto Null(u)
 Pc=I-P;
 % The vector A-orthogonal to v from W is obtained by projecting one of the
 % column vectors of W onto the space that is orth(u), in Null(u)
 projv=Pc*W(:,2);
 projv=projv/norm(projv);
 U(:,1:2)=[v,projv];

 for(k=3:Na);
    v=W(:,k);
    P=(A*U)*inv((A*U)'*(A*U))*(A*U)';
    Pc=(I-P);
    projv=Pc*W(:,k);
    U(:,k)=projv/norm(projv);
 end;