[U]=conjorth(W,A)

function [U]=conjorth(W,A)
%Compute a set of normalized, A-conjugate vectors from basis W.
%Vectors can be A-normalized outside of function easily.  See below.
%
%[U]=conjorth(W,A)
%
%INPUT
% W: A linear independet set of column vectors
% A: A full rank matrix.  An output is produced if A is not symmetric,
%    but A must be symmetric in order for U'*A*U=0;
%
%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. The columns of U are normalized.
%   
%   To A-normalize U(:,k), compute U(:,k)/(U(:,k)'*A*U(:,k))
%
%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);

% Make sure A and W are compatible sizes, and check for symmetry in A.
if(~isequal([Mw,Nw],[Ma,Na]))
    error('A and W must have identical dimensions');
end;

if(rank(A)<min([Ma,Na]))
    error('A must have full rank');
end;

if(rank(W)<min([Mw,Nw]))
    error('W must be full rank');
end;

if(norm(A-A') > 10^(-7))
    disp('A numerically asymmetric. U will be psuedo-A-conjugate only.')
    disp('This means transpose(U(:,2))*A*U(:,1) = 0, while:');
    disp('transpose(U(:,1))*A*U(:,2) ~= 0');
end;

% MAKING w1 and w2 A-conjugate to each other: Illustrating the method via
% construction of the first two conjugate-orthogonal vectors:
U=zeros(size(W));

% 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);

% Update U. Now v and projv are A-orthogonal so (A*v)'*projv = 0.
U(:,1:2)=[v,projv];

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