function [s,U,V,smlflag] = appsvd(a,b,tolsd)

% [s,U,V,smlflag] = appsvd(a,b,tolsd)

%

% SVD of a bidiagonal matrix

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

% using twisted factorization method.  Applied only when neccessary

% as small subdiagonals may be skipped over, resulting in a series of

% smaller SVDs.

%

% Inputs: 

%     a -- main diagonal of a bidiagonal matrix

%     b -- superdiagonal of a bidiagonal matrix

%     tolsd -- tolerance for small superdiagonal elements

%

% Outputs:

%     s -- singular (Takagi) value vector

%     U, V -- if no small diagonal is found than, 

%             left, right singular vector (Takagi vector) matrices

%              so that

%

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

%

%             otherwise the last U,V in a series of computed SVDs are 

%             returned.

%     smlflag -- set to 1 if small superdiagonal found, 0 otherwise

%

% Dependency:

%   gcstsvdt -- compute the svd of an upperbidiagonal matrix with

%               nonzero superdiagonal

%   clasv2 -- computes the svd of a 2x2 upper triangular matrix

%

% Kevin Browne McMaster Univ. June 2007

% Revised   April 2008, S. Qiao



n = length(a);



s = zeros(n, 1);

start_pos = 1;

smlflag = 0;



% 

k = 1;

while k <= n

    

    svdreq = 0; % SVD required?

    if (k == n) 

        svdreq = 1;

    elseif (abs(b(k)) < tolsd)

        fprintf('\nEncountered small superdiagonal at column %d', k+1);

        svdreq = 1;

        smlflag = 1; % a small superdiagonal has been found

    end

    

    if (svdreq == 1)



        % matrix order > 2

        if ((k - start_pos) > 2)

            

            % compute bidiagonal SVD, twisted factorization

            [s(start_pos:k),V,U] = gcstsvdt(a(start_pos:k),b(start_pos:k-1));

            

        % matrix order == 2x2

        elseif ((k - start_pos) > 0)

            

            % compute 2x2 upper triangular SVD            

            f = a(start_pos);

            g = b(start_pos);

            h = a(start_pos+1);

            [smin,smax,sr,cr,sl,cl] = clasv2(f,g,h);

            U = [sign(f),0;0,sign(h)] * [cl,-sl';sl,cl'];

            V = [cr,-sr';sr,cr'];



            % extract the singular values

            s(start_pos:k) = [smax,smin];

            

        end % if ((k - start_pos ...



        % Report matrix SVD position, size, error

        %B = diag(a(start_pos:k)) + diag(b(start_pos:k-1),1);

        %tmp_err = norm(B - (U * diag(s(start_pos:k)) * V'));

        %fprintf('\nComputed SVD of B %dx%d matrix',k+1-start_pos, ...

        %                                           k+1-start_pos);

        %fprintf('\nError in SVD: %d\n', tmp_err);

            

        k = k + 1;

        % skip over subsequent small superdiagonals

        while ((k <= (n-1)) & (abs(b(k)) < tolsd))

            s(k+1) = a(k+1);

            k = k + 1;

        end % while ((k <= ....

        

        % very last super diagonal was small? we are done

        if (k == n) 

            k = k + 1;

        end

        

        start_pos = k;   

        

    % keep looking for small subdiagonals

    else

        k = k + 1;

    end % if (svdreq == 1)

    

end % while k <= n



% be sure to order s properly

s = sort(s);

s(n:-1:1) = s(1:n);