%=======================================================================
% ch08_power_method.m
%
% Find the dominant eigenvalue and eigenvector by the power method.
%
% Origin:
%   MATLAB translation of ch08_power_method.f90 by Alexander Godunov,
%   based on Power.f90. The iteration y <- A*y / ||A*y|| is retained.
%   A signed Rayleigh quotient estimates the eigenvalue, and the
%   eigenpair residual tests convergence.
%
% Example:
%   A = [ 1  2  3; 2  2 -2; 3 -2  4 ], y0 = [ 1, 1, 1 ].
%   Its eigenvalues are approximately -2.54138127, 3.54138127,
%   and 6.00000000. The dominant eigenvalue in magnitude is 6.
%
% The method finds one eigenpair. It generally requires a unique
% dominant eigenvalue in magnitude and an initial vector with a
% nonzero component along its eigenvector. A small change in the
% eigenvalue estimate alone does not verify the eigenvector.
%=======================================================================

clear; clc;

A = [1.0,  2.0,  3.0;
     2.0,  2.0, -2.0;
     3.0, -2.0,  4.0];
y = [1.0; 1.0; 1.0];
eps = 1.0e-7;
max_iter = 1000;

fprintf('Power method for the dominant eigenpair\n');
fprintf('\nMatrix A:\n');
for i = 1:size(A,1)
    fprintf('%12.6f%12.6f%12.6f\n', A(i,:));
end
fprintf('\nInitial vector:%12.6f%12.6f%12.6f\n', y);
fprintf('Residual tolerance = %.4e\n', eps);

[y, lambda, iterations, residual, info] = ...
    power_method(A, y, eps, max_iter);
if info == 1
    error('Invalid input: check the initial vector and controls.');
elseif info == 2
    error('A*y is zero; the iteration cannot be normalized.');
end

fprintf('\nDominant eigenvalue = %12.6f\n', lambda);
fprintf('Unit eigenvector:\n');
for i = 1:numel(y)
    fprintf('y(%d) = %12.6f\n', i, y(i));
end
fprintf('Iterations = %d\n', iterations);
fprintf('Eigenpair residual ||A*y-lambda*y|| = %.4e\n', residual);
if info == 3
    error('Iteration limit reached before the tolerance was met.');
end


function [y, lambda, iterations, residual, info] = power_method(A, y, eps, max_iter)
%=======================================================================
% Approximate the eigenpair of largest eigenvalue magnitude.
%
% Input:
%   A        - square real matrix, left unchanged
%   y        - nonzero initial vector, returned as a unit vector
%   eps      - positive tolerance for the eigenpair residual
%   max_iter - maximum number of normalized power iterations
%
% Output:
%   lambda     - signed Rayleigh quotient y'*A*y / (y'*y)
%   iterations - completed power iterations
%   residual   - Euclidean norm of A*y - lambda*y
%   info       - 0 if residual <= eps; 1 for invalid input;
%                2 if A*y is zero; 3 at the iteration limit
%
% Normalize the starting vector. Each power step multiplies by A and
% normalizes. The Rayleigh quotient retains the sign of a negative
% dominant eigenvalue; ||A*y|| alone cannot do so. A small residual
% verifies the pair but cannot establish dominance if the initial
% vector misses the dominant eigendirection.
%=======================================================================

lambda = 0.0;
iterations = 0;
residual = Inf;
info = 1;
n = numel(y);
if n < 1 || ~isreal(A) || ~isreal(y) || ~isequal(size(A), [n,n]) ...
        || ~isscalar(eps) || eps <= 0 || ~isscalar(max_iter) ...
        || max_iter < 1 || max_iter ~= floor(max_iter)
    return
end
y = y(:);

% An eigenvector has no preferred length, so start with ||y|| = 1.
norm_y = norm(y);
if norm_y == 0.0
    return
end
y = y/norm_y;

for k = 1:max_iter
    % Power step: multiply by A, then normalize to prevent growth.
    z = A*y;
    norm_z = norm(z);
    if norm_z == 0.0
        info = 2;
        return
    end
    y = z/norm_z;

    % Estimate the signed eigenvalue for this normalized vector.
    Ay = A*y;
    lambda = y'*Ay;

    % Check the full eigenpair, not only the change in eigenvalue.
    residual = norm(Ay-lambda*y);
    iterations = k;
    if residual <= eps
        info = 0;
        return
    end
end
info = 3;
end