include iostream include iomanip include ios include vector include cm

#include <iostream>
#include <iomanip>
#include <ios>
#include <vector>
#include <cmath>
#include <functional>
#include <limits>
#include <sstream>

using namespace std;

typedef long double ld;
typedef vector<ld> VLD;
typedef vector<VLD> 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 MatrixPrint(const VVLD& a, const string& name = "", int precision = 3)
{
  cout << "==============\n";
  if (name != "") cout << name << '\n';
  for (const VLD& row : a)
  {
    cout << '[';
    for (const ld& x : row)
      cout << fixed << setprecision(precision) << x << ' ';
    cout << "]\n";
  }
  cout << "==============\n";
}

inline void VectorPrint(const VLD& a, const string& name = "", int precision = 3)
{
  cout << "==============\n";
  if (name != "") cout << name << '\n';
  cout << '[';
  for (const ld& x : a)
    cout << fixed << setprecision(precision) << x << ' ';
  cout << "]\n";
  cout << "==============\n";
}

inline void LuDecomposition(const VVLD& A, VVLD& L, VVLD& U)
{
  int n = A.size();
  L.assign(n, VLD(n, 0));
  U.assign(n, VLD(n, 0));
  for (int i = 0; i < n; ++i)
    U[i][i] = 1;

  for (int i = 0; i < n; ++i)
  {
    //i column of L
    for (int j = i; j < n; ++j) { //calculating L[j][i]
      L[j][i] = A[j][i];
      for (int t = 0; t < i; ++t)
        L[j][i] -= L[j][t] * U[t][i];
    }

    //i row of U
    for (int j = i + 1; j < n; ++j)
    {
      U[i][j] = A[i][j];
      for (int t = 0; t < i; ++t)
        U[i][j] -= L[i][t] * U[t][j];
      U[i][j] /= L[i][i];
    }
  }
}

inline void LuInverse(const VVLD& L, const VVLD& U, VVLD& X)
{
  int n = L.size();
  X.assign(n, VLD(n, 0));
  for (int i = n - 1; i >= 0; --i)
  {
    X[i][i] = 1;
    for (int k = i + 1; k < n; ++k)
      X[i][i] -= L[k][i] * X[i][k];
    X[i][i] /= L[i][i];

    for (int j = i - 1; j >= 0; --j)
      for (int k = j + 1; k < n; ++k)
        X[j][i] -= U[j][k] * X[k][i];

    for (int j = i - 1; j >= 0; --j){
      for (int k = j + 1; k < n; ++k)
        X[i][j] -= X[i][k] * L[k][j];
      X[i][j] /= L[j][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 CholDecomposition(const VVLD& A, VVLD& U)
{
  int n = A.size();
  U.assign(n, VLD(n, 0));
  for (int i = 0; i < n; ++i)
  {
    U[i][i] = A[i][i];
    for (int k = 0; k < i; ++k)
      U[i][i] -= U[k][i] * U[k][i];
    U[i][i] = sqrtl(U[i][i]);

    for (int j = i + 1; j < n; ++j)
    {
      U[i][j] = A[i][j];
      for (int k = 0; k < i; ++k)
        U[i][j] -= U[k][i] * U[k][j];
      U[i][j] /= U[i][i];
    }
  }
}

inline void CholSolve(const VVLD& U, const VLD& b, VLD& x) //U is result CholDecomposition(A)
{
  int n = U.size();
  x.assign(n, 0);
  for (int i = 0; i < n; ++i)
  {
    x[i] = b[i];
    for (int j = 0; j < i; ++j)
      x[i] -= U[j][i] * x[j];
    x[i] /= U[i][i];
  }
  for (int i = n - 1; i >= 0; --i)
  {
    for (int j = i + 1; j < n; ++j)
      x[i] -= U[i][j] * x[j];
    x[i] /= U[i][i];
  }
}

inline ld getNorm(const VLD& x)
{
  ld ans = 0;
  for (const ld& y : x)
    ans = max(ans, fabsl(y));
  return ans;
}
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<ld(const VLD&)> 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;
}

inline void SolveAll(const VVLD& A, const VLD& b, int outPrecision = 3)
{
  cout << "Ax=b" << endl;
  MatrixPrint(A, "Matrix A:", outPrecision);
  VectorPrint(b, "Vector B:", outPrecision);

  VVLD L, U;
  LuDecomposition(A, L, U);
  cout << "A = LU" << endl;
  MatrixPrint(L, "Matrix L:", outPrecision);
  MatrixPrint(U, "Matrix U:", outPrecision);

  VVLD invA;
  LuInverse(L, U, invA);
  MatrixPrint(invA, "Matrix A^(-1):", outPrecision);

  VLD x;
  MatrixMultiply(invA, b, x);
  VectorPrint(x, "Vector x = A^(-1)b (for eq. Ax=b):", outPrecision);

  cout << "Cholesky method. A = U'*U:" << endl;
  CholDecomposition(A, U);
  MatrixPrint(U, "Matrix U:", outPrecision);
  CholSolve(U, b, x);
  VectorPrint(x, "Vector x (for eq. U'Ux=b):", outPrecision);


  cout << "Pm method." << endl;
  ld eigValMax;
  VLD eigVecMax;
  PmAlgo(A, eigVecMax, eigValMax);
  stringstream ss;
  ss << "Max eigen value: " << fixed << setprecision(outPrecision) << eigValMax << "\nEigen vector: ";
  VectorPrint(eigVecMax, ss.str());
  ss.str("");
  ss.clear();

  VVLD tA = A;
  for (int i = 0; i < (int)A.size(); ++i)
    tA[i][i] -= eigValMax;

  ld eigValMin;
  VLD eigVecMin;
  PmAlgo(tA, eigVecMin, eigValMin);
  eigValMin += eigValMax;
  ss << "Min eigen value: " << fixed << setprecision(outPrecision) << eigValMin << "\nEigen vector: ";
  VectorPrint(eigVecMin, ss.str());
  ss.str().clear();
}

int main()
{
  /*VVLD A = { {1, 2, 3}, {2, 5, 6}, {3, 6, 11} };
  VLD b = {5, 6, 10};*/
  int n;
  cin >> n;
  VVLD A(n, VLD(n));
  VLD b(n);
  for (int i = 0; i < n; ++i)
    for (int j = 0; j < n; ++j)
      cin >> A[i][j];
  for (int i = 0; i < n; ++i)
    cin >> b[i];
  SolveAll(A, b);

  return 0;
}