function [iter,V,D] = alhqr(A,maxiter) % find the eigenvalues and corresponding eigenvectors of equation % (a-eval)*evect=0 using the power method with the shifted Hessenberg % QR algorithm. % based on a Matlab script by Gorka Erana Robles, June 27, 2017 % function [iter,V,D] = alhqr(A,maxiter) % (c) 2020 Alain Hebert, Polytechnique Montreal % Check input matrix size if ~isreal(A), error('Input matrix must be real'), end [m,n] = size(A); if (m ~= n), error('Input matrix must be square'), end % Initialize variables H = A; V = cell(1,n-2); % Perform Householder transformation to upper Hessenberg form for k = 1:(n-2) v = H((k+1):n,k); sgn = sign(v(1)); if (sgn == 0), sgn = 1; end v(1) = v(1) + sgn * norm(v); v = v / norm(v); H((k+1):n,k:n) = H((k+1):n,k:n) - 2 * v * (v' * H((k+1):n,k:n)); H(:,(k+1):n) = H(:,(k+1):n) - (2 * (H(:,(k+1):n) * v)) * v'; V{k} = v; end % Construct orthogonal matrix Q Q = eye(n); for j = (n-2):-1:1 Q((j+1):n,:) = Q((j+1):n,:) - (2 * V{j}) * (V{j}' * Q((j+1):n,:)); end % Perform Schur factorization i2 = n; c = zeros(1,n); s = zeros(1,n); iter = 0; while 1 iter = iter + 1; if (iter > maxiter) error('Maximum number of iterations exceeded.') end % Check subdiagonal for near zeros, deflating points. Finds deflating rows % on a complex Schur form matrix. i1 = i2; normH = norm(H,'fro'); while (i1 > 1) if (abs(H(i1,i1-1)) < eps*normH) H(i1,i1-1) = 0; if (i1 == i2) i2 = i1 - 1; i1 = i1 - 1; else break end else i1 = i1 - 1; end end % End the function if H is upper triangular if ( i2==1 ), break, end % Compute Wilkinson shift kappa = H(i2,i2); sum = abs(H(i2-1,i2-1)) + abs(H(i2-1,i2)) + abs(H(i2,i2-1)) + abs(H(i2,i2)); if sum ~= 0 q = (H(i2-1,i2)/sum)*(H(i2,i2-1)/sum); if q ~= 0 p = 0.5*((H(i2-1,i2-1)/sum) - (H(i2,i2)/sum)); r = sqrt(p*p + q); if ( (real(p)*real(r) + imag(p)*imag(r)) < 0 ), r = -r; end kappa = kappa - sum*(q/(p+r)); end end % Apply shift to the element of the diagonal that is left out of the loop H(i1,i1) = H(i1,i1) - kappa; for j = i1:i2-1 % Loop reducing the matrix to triangular form % Apply Givens rotation so that the subdiagonal is set to zero if H(j+1,j)==0 c(j) = 1; s(j) = 0; elseif H(j,j)==0 c(j) = 0; s(j) = 1; H(j,j) = H(j+1,j); H(j+1,j) = 0; else mu = H(j,j)/abs(H(j,j)); tau = abs(real(H(j,j))) + abs(imag(H(j,j))) + abs(real(H(j+1,j))) ... + abs(imag(H(j+1,j))); nu = tau*sqrt(abs(H(j,j)/tau)^2 + abs(H(j+1,j)/tau)^2); c(j) = abs(H(j,j))/nu; s(j) = mu*(H(j+1,j)')/nu; H(j,j) = nu*mu; H(j+1,j) = 0; end % Apply shift to diagonal H(j+1,j+1) = H(j+1,j+1) - kappa; % Modify the involved rows using a plane rotation t = c(j)*H(j,j+1:n) + s(j)*H(j+1,j+1:n); H(j+1,j+1:n) = c(j)*H(j+1,j+1:n) - (s(j)')*H(j,j+1:n); H(j,j+1:n) = t; end for k = i1:i2-1 % Loop applying the back multiplication using a plane rotation t = c(k)*H(1:k+1,k) + (s(k)')*H(1:k+1,k+1); H(1:k+1,k+1) = c(k)*H(1:k+1,k+1) - s(k)*H(1:k+1,k); H(1:k+1,k) = t; % Accumulate transformations using a plane rotation t = c(k)*Q(1:n,k) + (s(k)')*Q(1:n,k+1); Q(1:n,k+1) = c(k)*Q(1:n,k+1) - s(k)*Q(1:n,k); Q(1:n,k) = t; H(k,k) = H(k,k) + kappa; end H(i2,i2) = H(i2,i2) + kappa; end % Construct the orthonormal basis V=eye(n); for j=2:n for i=j-1:-1:1 denom=H(i,i)-H(j,j); if denom ~= 0 V(i,j)=V(i,j)-H(i,i+1:j)*V(i+1:j,j)/denom; end end end V=Q*V; D=diag(H); % Sort and normalize the eigensolution [Y,I_eig] = sort(abs(D)); I_eig = I_eig(n:-1:1); D = D(I_eig); V = V(:,I_eig); iset = find(abs(imag(D))>1.0e-10); index = iset(1:2:length(iset)); for i=1:n if ismember(i,index) j=iset(find(iset==i)+1); m=[V(i,i), V(i,j); V(j,i), V(j,j)]; m(:,1)=m(:,1)/norm(m(1:2,1)); m(:,2)=m(:,2)/norm(m(1:2,2)); AA=[real(m(1,1))+real(m(2,1)), -imag(m(1,1))-imag(m(2,1)); ... imag(m(1,1))+imag(m(2,1)), real(m(1,1))+real(m(2,1))]; BB=[real(m(1,2))+real(m(2,2)), -imag(m(1,2))-imag(m(2,2)); ... -imag(m(1,2))-imag(m(2,2)), -real(m(1,2))-real(m(2,2))]; CC=BB\AA; V(:,i)=V(:,i)*(CC(1,1)+sqrt(-1)*CC(1,2)); elseif ismember(i,iset) j=iset(find(iset==i)-1); if abs(D(i)-D(j)') > 1.0e-10, error('pathological ordering'), end D(i)=D(j)'; else D(i)=real(D(i)); end V(:,i)=V(:,i)/norm(V(:,i)); end D=diag(D);