function [d,V,U] = gcstsvdt(a,b)

% [d,V,U] = gcstsvdt(a,b)

%

% SVD of an upperbidiagonal matrix

%            B = diag(a) + diag(b,1)

% using twisted factorization method.

%

% Inputs: 

%     a -- main diagonal of an upperbidiagonal matrix

%     b -- superdiagonal of an upperbidiagonal matrix

%

% Outputs:

%     d -- singular (Takagi) value vector

%  U, V -- left, right singular vector (Takagi vector) matrices

%          so that

%              diag(a) + diag(b,1) = U*diag(d)*V.'

%

% Dependency:

%     gtwist -- compute the corresponding eigenvector using the 

%               twisted factorization.

%

% Reference:

%  W. Xu and S. Qiao.

%  A divide-and-conquer method for the Takagi factorization.

%  Technical Report No. CAS 05-01-SQ. Department of Computing and

%  Software, McMaster University, Hamilton, Ontario L8S 4K1, Canada.

%  February 2005.

%

% W. Xu and S. Qiao  McMaster Univ. May 2007

% Revised by Kevin Browne McMaster Univ. June 2007

% 



% get dimension of the matrix Q

n = length(a);



% get the singular values of a complex upperbidiagonal B, squared to get

% eigenvalues of B'*B

e = svd((diag(a) + diag(b,1))).^2;



% get the corresponding eigenvector for the eigenvalue, LQ decomposition

for k = 1:n

    V(:,k) = gtwist(a,b,e(k));

end



% convert eigenvalues of B'*B into singular values of B 

d = sqrt(e);



% effeciently compute U = B * V

for k=1:n

    for j=1:n-1

        U(j,k) =  V(j,k) * a(j) + V(j+1,k) * b(j);

    end

    U(n,k) = V(n,k) * a(n);  

    

    % normalize the columns of B * V

    U(:,k) = U(:,k) / norm(U(:,k));

end