typedef long double ld; typedef vector VLD; typedef vector VVLD; inline void MatrixMultiply(const VVLD& A, const VVLD& B, VVLD& result) { int a = A.size(); int b = A[0].size(); if (b != (int)B.size()) return; int c = B[0].size(); result.assign(a, VLD(c, 0)); for (int i = 0; i < a; ++i) for (int j = 0; j < c; ++j) for (int k = 0; k < b; ++k) result[i][j] += A[i][k] * B[k][j]; } inline void MatrixMultiply(const VVLD& A, const VLD& b, VLD& res) { int n = A.size(); res.assign(n, 0); for (int i = 0; i < n; ++i) for (int j = 0; j < n; ++j) res[i] += A[i][j] * b[j]; } inline void Mult(const VLD& a, ld x, VLD& res) { res.resize(a.size()); for (size_t i = 0; i < a.size(); ++i) res[i] = a[i] * x; } inline void Divide(const VLD& a, ld x, VLD& res) { res.resize(a.size()); for (size_t i = 0; i < a.size(); ++i) res[i] = a[i] / x; } inline void PmAlgo( const VVLD& A, VLD& eigVec, ld& eigVal, const function norm = getNorm, const ld EPS = 1e-5, const int maxIters = 10000) { int n = A.size(); VLD y(n, 1), yk(n); VLD x(n), xk(n); VLD lambda(n, 1e100), lambdak(n); Divide(y, norm(y), x); bool finded = false; for (int iter = 0; iter < maxIters && !finded; ++iter) { MatrixMultiply(A, x, yk); Divide(yk, norm(yk), xk); ld maxdelta = 0; for (int i = 0; i < n; ++i){ lambdak[i] = yk[i] / x[i]; maxdelta = max(maxdelta, fabsl(lambda[i] - lambdak[i])); } finded = (maxdelta <= EPS); swap(x, xk); swap(y, yk); swap(lambda, lambdak); } eigVal = 0; for (ld t : lambda) eigVal += t; eigVal /= n; eigVec = x; }