function [a,b,V,U] = gLanMPO(A,rv) % [a,b,V,U] = gLanMPO(A,rv) % % General Lanczos bidiagonalization with modified partial orthogonalization. % Given a (complex) matrix A, this computes % a, b, and unitary V and U such that % A*V = U*(diag(a) + diag(b,1)) % % Input % A m-by-n (m >= n) complex matrix % 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 % % Reference % Gene H. Golub and Charles F. Van Loan. % Matrix Computations, 3rd Edition. % The Johns Hopkins University Press, % Baltimore, MD. 1996. p. 495. [m,n] = size(A); IM = sqrt(-1); % constants SQRTEPS = sqrt(eps); MIDEPS = sqrt(SQRTEPS)^3; % between sqrt(eps) and eps a = zeros(n,1); b = zeros(n-1,1); V = zeros(n,n); U = zeros(m,n); p = rv/norm(rv); % initialize a unit 2-norm vector bOld = 1.0; uOld = zeros(m,1); 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; % no of vectors selected for ortho upV(1) = 0; lowV(1) = 0; upU(1) = 0; lowU(1) = 0; interNumV = 0; interNumU = 0; for j=1:n vj = p/bOld; V(:,j) = vj; % find a column of U r = A*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-1) if (abs(muCur(k)) > SQRTEPS) % lost orthogonality interNumU = interNumU + 1; % find the upper bound curpos = k + 1; while ((curpos <= (j - 1)) & (abs(muCur(curpos)) >= MIDEPS)) curpos = curpos + 1; % nearly lost orthogonality end % while upU(interNumU) = curpos - 1; % find the lower bound curpos = k - 1; while ((curpos > 0) & (abs(muCur(curpos)) >= MIDEPS)) curpos = curpos - 1; % nearly lost orthogonality 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 orthogonalization 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 = A'*uj - a(j)*vj; b(j) = norm(p); % 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; % not the second time, determine intervals if (secondV == 0) interNumV = 0; k = 1; while k <= j if (abs(nuCur(k)) > SQRTEPS) % lost orthogonality interNumV = interNumV + 1; % find the upper bound curpos = k + 1; while ((curpos <= j) & (abs(nuCur(curpos)) >= MIDEPS)) curpos = curpos + 1; % nearly lost orthogonality end % while upV(interNumV) = curpos - 1; % find the lower bound curpos = k - 1; while ((curpos > 0) & (abs(nuCur(curpos)) >= MIDEPS)) curpos = curpos - 1; % nearly lost orthogonality 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 % now we have intervals, carry out orthogonalization if (interNumV > 0 | (secondV == 1)) for (k = 1:interNumV) % for each interval for (i = lowV(k):upV(k)) % do orthogonalization p = p - (V(:,i)'*p)*V(:,i); nuCur(i) = eps*1.5*(randn + IM*randn); % reset ortho estimates end % for i nVecV = nVecV + upV(k) - lowV(k) + 1; % count the vectors selected if (secondV == 1) % this is the second time lowV(k) = 0; upV(k) = 0; else % this is the first time % adjust orthogonalization intervals for the second time lowV(k) = max(1, lowV(k) - 1); upV(k) = min(j, upV(k) + 1); end % if second == 1 end % for k secondV = ~secondV; % repeat if not already second time b(j) = norm(p); % recalculate b(j) end % if interNum > 0 ... bOld = b(j); end % if j