%
% Matlab code prepared in August 2003 by Brittainy Owens
%   and Gary Howell at ERDC to send to Jesse Barlow 
%   in answer to his request for a demonstration of how
%   blocking could be applied to the
%   Ralha-Barlow one sided orthogonalization algorithm for
%   bidiagonalization. 
%   
%   Brittainy was then a summer intern at ERDC and a sophomore
%     at Hampton Roads.
%
%   Comments added by GWH April 13, 2006. 
%     Running the code also requires house.m 
%  
  
b=eye(n)
xold=x
g(1)=norm(x(:,1));
u(:,1)=x(:,1)/g(1);
z1=x(:,2:n)'*u(:,1);
temp=z1;
house;
x(:,2:n)= x(:,2:n)-x(:,2:n)*temp*temp'; 

%  This is for the special case that blocks are of size 2.

for k=2:n-2:2
  p(k)=u(:,k-1)'*x(:,k);
  x(:,k)=x(:,k)-p(k)*u(:,k-1);
  g(k)=norm(x(:,k));u(:,k)=x(:,k)/g(k);
  zk=x(:,k+1:n)'*u(:,k);
  temp=zk;
  house;
  tempold=temp

  x(:,k+1)=x(:,k+1)-x(:,k+1:n)*(temp*temp(1))
  p(k+1)=u(:,k)'*x(:,k+1);
  x(:,k+1)=x(:,k+1)-p(k+1)*u(:,k);
  g(k+1)=norm(x(:,k+1));u(:,k+1)=x(:,k+1)/g(k+1);
  zk=x(:,k+2:n)'*u(:,k+1);
  temp=zk;
  house;  	
  y=[tempold(2:n-k) temp]
  x(:,k+2:n)=x(:,k+2:n)-x(:,k+2:n)*y*y';

end

if(n>2)
  p(n-1)=u(:,n-2)'*x(:,n-1);
  x(:,n-1)=x(:,n-1)-p(n-1)*u(:,n-2); 
  g(n-1)=norm(x(:,n-1));u(:,n-1)=x(:,n-1)/g(n-1);
end 
if(n>1)
  p(n)=u(:,n-1)'*x(:,n);
  x(:,n)=x(:,n)-p(n)*u(:,n-1);  
  g(n)=norm(x(:,n));u(:,n)=x(:,n)/g(n);
end

for i=1:n
  b(i,i)=g(i);
end
for i=1:n-1
  b(i,i+1)=p(i+1);
end
err=svd(b)-svd(xold)