function [iter,V,D] = almeig(A,m,eps) % find the first m eigenvalues and corresponding eigenvectors of % equation (a-eval)*evect=0 using the power method with the method % of orthogonal iteration. % function [iter,V,D] = almeig(A,m,eps) % (c) 2020 Alain Hebert, Polytechnique Montreal n=size(A,1); if m > n, error('almeig: m>n forbidden.'), end Q=tril(ones(n,m)); iter=0; while true iter=iter+1; if iter > 500, error('almeig: unable to converge.'), end GAR=Q; Z=A*Q; [X, tau] = alst2c(Z); % economy-size QR factorization R = triu(X(1:m,1:m)); Q = eye(n,m); for j = m:-1:1 w = [1; X(j+1:n,j)]; Q(j:n,:) = Q(j:n,:)+tau(j)*w*(w'*Q(j:n,:)); end err1=max(abs(Q),[],'all'); err2=max(abs(GAR-Q),[],'all'); if err2 <= err1*eps, break, end end % construct the orthonormal basis V=eye(n,m); for j=2:m for i=j-1:-1:1 denom=R(i,i)-R(j,j); if denom ~= 0 V(i,j)=V(i,j)-R(i,i+1:j)*V(i+1:j,j)/denom; end end end V=Q*V(1:m,:); D=diag(diag(R(1:m,1:m)));