/*=======================================================================
  ch08_power_method.cpp

  Find the dominant eigenvalue and eigenvector by the power method.

  Origin:
    C++ 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.
=======================================================================*/

#include <cmath>
#include <iomanip>
#include <iostream>
#include <vector>

using Matrix = std::vector<std::vector<double>>;

void matrix_vector(const Matrix& a, const std::vector<double>& v,
                   std::vector<double>& w);
int power_method(const Matrix& a, std::vector<double>& y, double eps,
                 int max_iter, double& lambda, int& iterations,
                 double& residual);

int main()
{
    const Matrix a = {{1.0,  2.0,  3.0},
                      {2.0,  2.0, -2.0},
                      {3.0, -2.0,  4.0}};
    std::vector<double> y = {1.0, 1.0, 1.0};
    const double eps = 1.0e-7;
    const int max_iter = 1000;

    std::cout << "Power method for the dominant eigenpair\n"
              << "\nMatrix A:\n" << std::fixed << std::setprecision(6);
    for (const auto& row : a) {
        for (double value : row)
            std::cout << std::setw(12) << value;
        std::cout << '\n';
    }
    std::cout << "\nInitial vector:";
    for (double value : y)
        std::cout << std::setw(12) << value;
    std::cout << "\nResidual tolerance = " << std::scientific
              << std::setprecision(4) << eps << '\n';

    double lambda, residual;
    int iterations;
    const int info = power_method(a, y, eps, max_iter,
                                  lambda, iterations, residual);
    if (info == 1) {
        std::cerr << "\nInvalid input: check the initial vector and controls.\n";
        return 1;
    }
    if (info == 2) {
        std::cerr << "\nA*y is zero; the iteration cannot be normalized.\n";
        return 1;
    }

    std::cout << "\nDominant eigenvalue = " << std::fixed
              << std::setprecision(6) << std::setw(12) << lambda << '\n'
              << "Unit eigenvector:\n";
    for (std::size_t i = 0; i < y.size(); ++i)
        std::cout << "y(" << i + 1 << ") = " << std::setw(12) << y[i] << '\n';
    std::cout << "Iterations = " << iterations << '\n'
              << "Eigenpair residual ||A*y-lambda*y|| = "
              << std::scientific << std::setprecision(4) << residual << '\n';
    if (info == 3) {
        std::cerr << "Iteration limit reached before the tolerance was met.\n";
        return 1;
    }
    return 0;
}

int power_method(const Matrix& a, std::vector<double>& y, double eps,
                 int max_iter, double& lambda, int& iterations,
                 double& residual)
/*=======================================================================
  Approximate the eigenpair of largest eigenvalue magnitude.

  Input:
    a        - square real matrix, left unchanged
    y        - nonzero initial vector, replaced by a unit vector
    eps      - positive tolerance for the eigenpair residual
    max_iter - maximum number of normalized power iterations

  Output:
    lambda     - signed Rayleigh quotient y^T*A*y / (y^T*y)
    iterations - completed power iterations
    residual   - Euclidean norm of A*y - lambda*y
  Return value:
    0 if residual <= eps; 1 for invalid input; 2 if A*y is zero;
    3 if the iteration limit is reached.

  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 returned eigenpair, but cannot establish dominance if
  the initial vector misses the dominant eigendirection.
=======================================================================*/
{
    lambda = 0.0;
    residual = 0.0;
    iterations = 0;
    const std::size_t n = y.size();
    if (n == 0 || a.size() != n || eps <= 0.0 || max_iter < 1)
        return 1;
    for (const auto& row : a)
        if (row.size() != n)
            return 1;

    double norm_y = 0.0;
    for (double value : y)
        norm_y += value * value;
    norm_y = std::sqrt(norm_y);
    if (norm_y == 0.0)
        return 1;
    for (double& value : y)
        value /= norm_y;

    std::vector<double> z(n), ay(n);
    for (int k = 1; k <= max_iter; ++k) {
        // Power step: multiply by A, then normalize to prevent growth.
        matrix_vector(a, y, z);
        double norm_z = 0.0;
        for (double value : z)
            norm_z += value * value;
        norm_z = std::sqrt(norm_z);
        if (norm_z == 0.0)
            return 2;
        for (std::size_t i = 0; i < n; ++i)
            y[i] = z[i] / norm_z;

        // Estimate the signed eigenvalue for this normalized vector.
        matrix_vector(a, y, ay);
        lambda = 0.0;
        for (std::size_t i = 0; i < n; ++i)
            lambda += y[i] * ay[i];

        // Check the full eigenpair, not only the change in eigenvalue.
        double residual_sq = 0.0;
        for (std::size_t i = 0; i < n; ++i)
            residual_sq += (ay[i] - lambda*y[i]) * (ay[i] - lambda*y[i]);
        residual = std::sqrt(residual_sq);
        iterations = k;
        if (residual <= eps)
            return 0;
    }
    return 3;
}

void matrix_vector(const Matrix& a, const std::vector<double>& v,
                   std::vector<double>& w)
/* Compute w=A*v explicitly, without calling a matrix library. */
{
    for (std::size_t i = 0; i < v.size(); ++i) {
        w[i] = 0.0;
        for (std::size_t j = 0; j < v.size(); ++j)
            w[i] += a[i][j] * v[j];
    }
}