function [a,b,V,U] = gFHLanMPOR(col,row,rv)
% [a,b,V,U] = gFHLanMPOR(col,row,rv)
%
% General Fast Hankel Lanczos bidiagonalization with modified
% partial orthogonalization and restart. Given a (complex) Hankel
% matrix A=hankel(col,row), this computes a, b, and unitary V and U
% such that
% A*V = U*(diag(a) + diag(b,1))
%
% Input
% col first column of Hankel
% row last row of Hankel
% rv n-by-1, start vector for V
% Outputs
% a,b diagonal and subdiagonal
% V n-by-n unitary matrix
% U m-by-n with orthnormal columns
% S. Qiao McMaster Univ. May 2007
% Revised by Kevin Browne McMaster Univ. June 2007
% Nov 2007
%
% Reference
% Gene H. Golub and Charles F. Van Loan.
% Matrix Computations, 3rd Edition.
% The Johns Hopkins University Press,
% Baltimore, MD. 1996. p. 495.
%
% Dependency
% ./fhmvmtconv.m Hankel to Toeplitz-ciruclant conversion
% ./fhmvmmultiply.m Fast general Hankel matrix-vector multiplication
% ./hankelfnrm.m Frobenius norm of an m-by-n (m>=n) Hankel matrix
m = length(col);
n = length(row);
% is A complex?
if ((imag(col) ~= 0) | (imag(row) ~= 0))
isComplex = 1;
else
isComplex = 0;
end
% convert Hankel to Toeplitz-circulant
% to be used in fast Hankel matrix-vector multiplication
cc = fhmvmtconv(col,row);
% convert complex-conjugate and transpose of Hankel to Toeplitz-circulant
% to be used in fast Hankel matrix-vector multiplication
colt = conj(col(1:n));
rowt = [conj(col(n:(m-1))); conj(row(1:n))];
cct = fhmvmtconv(colt,rowt);
% constants
IM = sqrt(-1);
SQRTEPS = sqrt(eps);
MIDEPS = sqrt(SQRTEPS)^3; % between sqrt(eps) and eps
RSTOL = sqrt(eps)*(hankelfnrm(col,row)/(m*n)); % tolerance for restart
MAXRS = 50; % maximal number of restarts
a = zeros(n,1);
b = zeros(n-1,1);
V = zeros(n,n);
U = zeros(m,n);
if (norm(rv) == 0)
rv = ones(n,1) - 2*rand(n,1);
if (isComplex == 1)
rv = rv + IM*(ones(n,1) - 2*rand(n,1));
end
end
p = rv/norm(rv); % initialize a unit 2-norm vector
bOld = 1.0; uOld = zeros(m,1);
normp = 1.0;
nuOld = zeros(n,1); % initial orthogonality estimates
nuOld(1) = 1.0;
nuCur = zeros(n,1);
muOld = zeros(n,1);
muCur = zeros(n,1);
muCur(1) = eps*m*0.6*(randn + IM*randn);
muCur(2) = 1.0;
secondU = 0; secondV = 0; % if do ortho in the next iteration
nVecV = 0; nVecU = 0; % number of vectors selected for orthog
upV(1) = 0; lowV(1) = 0; upU(1) = 0; lowU(1) = 0;
%
for j=1:n
vj = p / normp;
V(:,j) = vj;
% find a column of U
r = fhmvmultiply(cc,vj) - bOld*uOld;
a(j) = norm(r);
% estimate the orthogonality of r (U(:,j)) against U(:,1:j-1)
if (j > 2)
muOld(1:j-2) = (b(1:j-2).*nuCur(2:j-1) ...
+ a(1:j-2).*nuCur(1:j-2) ...
- b(j-1)*muCur(1:j-2))/a(j) ...
+ eps*(b(1:j-2) ...
+ a(j)*ones(j-2,1))*0.3*(randn + IM*randn);
% swap muOld(1:j-2) and muCur(1:j-2)
tmp = muOld(1:j-2);
muOld(1:j-2) = muCur(1:j-2);
muCur(1:j-2) = tmp;
muOld(j-1) = 1.0;
muCur(j-1) = eps*m*((a(2)*b(1))/(a(j)*b(j-1))) ...
*0.6*(randn + IM*randn);
muCur(j) = 1.0;
end % if (jc > 2)
% not the second time, determine intervals
if (secondU == 0)
interNumU = 0;
k = 1;
while k <= j
if (abs(muCur(k)) > SQRTEPS) % lost orthogonality
interNumU = interNumU + 1;
% find the upper bound
curpos = k + 1;
% nearly lost orthogonality
while ((curpos <= (j - 1)) & (abs(muCur(curpos)) >= MIDEPS))
curpos = curpos + 1;
end % while
upU(interNumU) = curpos - 1;
% find the lower bound
curpos = k - 1;
% nearly lost orthogonality
while ((curpos > 0) & (abs(muCur(curpos)) >= MIDEPS))
curpos = curpos - 1;
end % while
lowU(interNumU) = curpos + 1;
%
k = upU(interNumU) + 1; % continue search
else
k = k + 1;
end % if lost orthogonality
end % while k
end % if not second time
% now we have intervals, carry out orthogonalization
if (interNumU > 0 | (secondU == 1))
for (k = 1:interNumU) % for each interval
for (i = lowU(k):upU(k)) % do orthogonalization
r = r - (U(:,i)'*r)*U(:,i);
muCur(i) = eps*1.5*(randn + IM*randn); % reset ortho estimates
end % for i
nVecU = nVecU + upU(k) - lowU(k) + 1; % count the vectors selected
if (secondU == 1) % this is the second time
lowU(k) = 0; upU(k) = 0;
else % this is the first time
% adjust orthog intervals for the second time
lowU(k) = max(1, lowU(k) - 1);
upU(k) = min(j, upU(k) + 1);
end % if second == 1
end % for k
secondU = ~secondU; % repeat if not already second time
a(j) = norm(r); % recalculate a(j)
end % if interNum > 0 ...
uj = r/a(j);
uOld = uj;
U(:,j) = uj;
% find a column of V
if (j < n)
p = fhmvmultiply(cct,uj) - a(j)*vj;
normp = norm(p);
b(j) = normp;
% determine orthogonalization intervals
nRS = 0;
nEntry = 0; % number of entries into the following loop
while ((nEntry < 1) | ((normp < RSTOL) & (nRS < MAXRS)))
nEntry = nEntry + 1;
if (normp < (RSTOL)) % small subdiagonal detected
% orthogonalize p against all previous V(:,1:j)
interNumV = 1; % one ortho interval [1:j]
lowV(interNumV) = 1;
upV(interNumV) = j;
p = ones(n,1) - 2*rand(n,1); % restart
if (isComplex == 1)
p = p + IM*(ones(n,1) - 2*rand(n,1));
end
nRS = nRS + 1;
% update nuOld and nuCur
if (j >=2)
nuOld(1:j-1) = nuCur(1:j-1);
nuOld(j) = 1.0;
end
nuCur(j+1) = 1.0;
else % subdiagonal is not small, no restart
% estimate the orthogonality of p (V(:,j+1)) against V(:,1:j)
if (j > 2)
nuOld(1) = (a(1)'*muCur(1) - a(j)*nuCur(1))/b(j) ...
+ eps*(a(1) + b(j))*0.3*(randn + IM*randn);
nuOld(2:j-1) = (a(2:j-1).*muCur(2:j-1) ...
+ b(1:j-2).*muCur(1:j-2) ...
- a(j)*nuCur(2:j-1))/b(j) ...
+ eps*(a(2:j-1) ...
+ b(j)*ones(j-2,1))*0.3*(randn + IM*randn);
% swap nuOld(1:j-1) and nuCur(1:j-1)
tmp = nuOld(1:j-1);
nuOld(1:j-1) = nuCur(1:j-1);
nuCur(1:j-1) = tmp;
nuOld(j) = 1.0;
end % if j>2
nuCur(j) = eps*n*((a(1)*b(1)')/(a(j)*b(j)')) ...
*0.6*(randn + IM*randn);
nuCur(j+1) = 1.0;
if (secondV == 0) % not the second time, determine intervals
interNumV = 0;
k = 1;
while k <= (j-1)
if (abs(nuCur(k)) > SQRTEPS) % lost orthogonality
interNumV = interNumV + 1;
% find the upper bound
curpos = k + 1;
% nearly lost orthogonality
while ((curpos < (j + 1)) & ...
(abs(nuCur(curpos)) >= MIDEPS))
curpos = curpos + 1;
end % while
upV(interNumV) = curpos - 1;
% find the lower bound
curpos = k - 1;
% nearly lost orthogonality
while ((curpos > 0) & ...
(abs(nuCur(curpos)) >= MIDEPS))
curpos = curpos - 1;
end % while
lowV(interNumV) = curpos + 1;
%
k = upV(interNumV) + 1; % continue search
else
k = k + 1;
end % if lost orthogonality
end % while k
end % if not second time
end
% now we have intervals, carry out orthogonalization
if (interNumV > 0 | (secondV == 1))
% for each interval, do orthog, reset ortho estimates
for (k = 1:interNumV)
for (i = lowV(k):upV(k))
p = p - (V(:,i)'*p)*V(:,i); % do orthogonalization
nuCur(i) = eps*1.5*(randn + IM*randn);
end % for i
nVecV = nVecV + upV(k) - lowV(k) + 1;
% count the number of vectors selected
if (secondV == 1) % this is the second time
lowV(k) = 0; upV(k) = 0;
else
% adjust orthog intervals for the second time
lowV(k) = max(1, lowV(k) - 1);
upV(k) = min(j + 1, upV(k) + 1);
end % if
end % for k
secondV = ~secondV; % repeat if not already second time
normp = norm(p);
end % if
end
bOld = b(j);
if (normp < RSTOL)
fprintf('\nRestart %d times', nRSV);
a = a(1:j);
b = b(1:j-1);
nsteps = j;
return
end
% if no restart, we may have performed orthogonalization,
% recalculate b(j); if restart, we use the old small b(j)
% instead of normr.
if (nRS == 0) % no restart
b(j) = normp; % recalculate b(j)
end
end % if j<n
end % for j
%fprintf('\n %d U vectors selected for orthogonalization.', nVecU)
%fprintf('\n %d V vectors selected for orthogonalization.', nVecV)