for_gui_mmn_do

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

extern "C" {
double f_test(double x, double y) {
  double T = M_PI * x * y;
  double sin2T = std::pow(std::sin(T), 2);
  double cos2T = std::pow(std::cos(T), 2);
  double y2_x2 = y * y + x * x;
  return -2 * M_PI * std::exp(std::pow(std::sin(T), 2)) *
         (2 * M_PI * y2_x2 * cos2T * sin2T - M_PI * y2_x2 * sin2T +
          M_PI * y2_x2 * cos2T);
}
}

extern "C" {
double f_main(double x, double y) {
  return -std::pow(std::sin(M_PI * y * x), 2);
}
}
extern "C" {
double u(double x, double y) {
  return std::exp(std::pow(std::sin(M_PI * x * y), 2));
}
}
extern "C" {
void Write_test(std::vector<std::vector<double>> v, const int n, const int m,
                const int S, const double eps_max, const double r_max,
                const double z_max) // функция записи массива в бинарный файл
{
  std::ofstream in("/home/syatov430/VAZHNO/NM_labs_3kurs_2semestr/NM_lab2/lab2/"
                   "data_test.txt"); // подключение текстового файла для записи
                                     // в бинарном режиме
  for (int i = 0; i < n; i++) {
    for (int j = 0; j < m; j++) {
      if (j != m - 1)
        in << v[i][j] << ", "; // запись в файл поэлементно
      else
        in << v[i][j] << "\n";
    }
  }
  in.close(); // закрываем файл

  std::ofstream in1("/home/syatov430/VAZHNO/NM_labs_3kurs_2semestr/NM_lab2/"
                    "lab2/z_it.txt"); // подключение текстового файла для записи
                                      // в бинарном режиме
  for (int i = 0; i < n; i++) {
    for (int j = 0; j < m; j++) {
      if (j != m - 1)
        in1 << fabs(v[i][j] - u(i * (1 / (n - 1)), j * (1 / (m - 1))))
            << ", "; // запись в файл поэлементно
      else
        in1 << fabs(v[i][j] - u(i * (1 / (n - 1)), j * (1 / (m - 1)))) << "\n";
    }
  }
  in1.close(); // закрываем файл

  std::ofstream in2("/home/syatov430/VAZHNO/NM_labs_3kurs_2semestr/NM_lab2/"
                    "lab2/param_test.txt"); // подключение текстового файла для
                                            // записи в бинарном режиме
  in2 << S << ", " << eps_max << ", " << r_max << ", " << z_max;
  in2.close(); // закрываем файл
}
}

extern "C" {
void Write_main_v(std::vector<std::vector<double>> v, const int n, const int m,
                  const double v_v2, const int n2,
                  const int m2) // функция записи массива в бинарный файл
{
  std::ofstream in("/home/syatov430/VAZHNO/NM_labs_3kurs_2semestr/NM_lab2/lab2/"
                   "data_main.txt"); // подключение текстового файла для записи
                                     // в бинарном режиме
  for (int i = 0; i < n; i++) {
    for (int j = 0; j < m; j++) {
      if (j != m - 1)
        in << v[i][j] << ", "; // запись в файл поэлементно
      else
        in << v[i][j] << "\n";
    }
  }
  in.close(); // закрываем файл
  std::ofstream in2(
      "/home/syatov430/VAZHNO/NM_labs_3kurs_2semestr/NM_lab2/lab2/"
      "param_main_v2.txt"); // подключение текстового файла для записи в
                            // бинарном режиме
  in2 << v_v2 << n2 << m2;
  in2.close(); // закрываем файл
}
}

void Write_main_param(
    const int n, const int m, const int S, const double eps_max,
    const double r_max) // функция записи массива в бинарный файл
{

  std::ofstream in2("/home/syatov430/VAZHNO/NM_labs_3kurs_2semestr/NM_lab2/"
                    "lab2/param_main.txt"); // подключение текстового файла для
                                            // записи в бинарном режиме
  in2 << n << m << S << ", " << eps_max << ", " << r_max;
  in2.close(); // закрываем файл
}
extern "C" {
int MinimalResiduals_test(const int nmax = 1000, const double _eps = 0.0000001,
                          const int _n = 3, const int _m = 3) {
  int Nmax = nmax; // максимальное число итераций (не менее 1)
  std::vector<std::vector<double>> v; // сеточная функция ѵ (x, y)
  v.resize(_n + 1);
  for (size_t i = 0; i < _n + 1; i++) {
    v[i].resize(_m + 1);
    for (size_t j = 0; j < _m + 1; j++) {
      v[i][j] = 0;
    }
  }
  std::vector<std::vector<double>> r_vec; // вектор невязки
  r_vec.resize(_n + 1);
  for (size_t i = 0; i < _n + 1; i++) {
    r_vec[i].resize(_m + 1);
    for (size_t j = 0; j < _m + 1; j++) {
      r_vec[i][j] = 0;
    }
  }
  std::vector<std::vector<double>> ar; // вектор невязки
  ar.resize(_n + 1);
  for (size_t i = 0; i < _n + 1; i++) {
    ar[i].resize(_m + 1);
    for (size_t j = 0; j < _m + 1; j++) {
      ar[i][j] = 0;
    }
  }
  double r1 = 0;
  int S = 0;         // счетчик итераций
  double eps = _eps; // заданная точность
  double eps_max = 0; // точность, достигнутая на текущей итерации
  double eps_cur = 0; // для подсчета точности на текущей итерации
  double a2, k2, h2; // ненулевые элементы матрицы (-А)
  const int n = _n, m = _m; // размерность сетки
  double r_max = 0;
  double a = 0;
  double b = 1;
  double c = 0;
  double d = 1; // границы области определения уравнения
  int i, j; // индексы
  double r = 0;
  double v_old; // старое значение преобразуемой компоненты вектора
  double v_new; // новое значение преобразуемой компоненты вектора и
  bool flag = false; // условие остановки
  double h = ((b - a) / n);
  double k = ((d - c) / m);
  double z = 0;
  double z_max = 0;
  double tau_s = 0;
  double num = 0; // tau=num/denominator
  double denominator = 1;
  h2 = -std::pow((n / (b - a)), 2);
  k2 = -std::pow((m / (d - c)), 2);
  a2 = -2 * (h2 + k2);
  for (int i = 0; i < n + 1; i++) {
    v[i][0] = u(a + i * h, a);
    v[i][m] = u(a + i * h, b);
  }
  for (int j = 0; j < m + 1; j++) {
    v[0][j] = u(c, c + j * k);
    v[n][j] = u(d, c + j * k);
  }
  do {
    z = 0;
    z_max = 0;
    eps_max = 0;
    r_max = -20;
    num = 0;
    denominator = 0;
    for (j = 1; j < m; j++) {
      for (i = n - 1; i > 0; i--) {
        r_vec[i][j] = f_test(a + i * h, c + j * k) -
                      (h2 * (v[i + 1][j] + v[i - 1][j]) +
                       k2 * (v[i][j + 1] + v[i][j - 1]) + v[i][j] * a2);
      }
    }
    for (j = 1; j < m; j++) {
      for (i = n - 1; i > 0; i--) {
        ar[i][j] =
            (h2 * (r_vec[i + 1][j] + r_vec[i - 1][j]) +
             k2 * (r_vec[i][j + 1] + r_vec[i][j - 1]) + r_vec[i][j] * a2);
        num += ar[i][j] * r_vec[i][j];
        denominator += ar[i][j] * ar[i][j];
      }
    }
    tau_s = num / denominator;
    for (j = 1; j < m; j++) {
      for (i = n - 1; i > 0; i--) {
        v_old = v[i][j];
        v_new = v_old + tau_s * r_vec[i][j];
        eps_cur = std::fabs(v_old - v_new);
        if (eps_cur > eps_max) {
          eps_max = eps_cur;
        }
        v[i][j] = v_new;
        r = fabs(f_test(a + i * h, c + j * k) -
                 (h2 * (v[i + 1][j] + v[i - 1][j]) +
                  k2 * (v[i][j + 1] + v[i][j - 1]) + v_old * a2));
        if (r > r_max)
          r_max = r;
        z = u(a + i * h, c + j * k) - v[i][j];
        if (z > z_max)
          z_max = z;
      }
    }

    S = S + 1;
    if ((eps_max <= eps) or (S >= Nmax)) {
      flag = true;
    }
  } while (!flag);
  Write_test(v, n + 1, m + 1, S, eps_max, r_max, z_max);
  return 0;
}
}

extern "C" {
std::vector<std::vector<double>>
MinimalResiduals_main(const int nmax = 1000, const double _eps = 0.0000001,
                      const int _n = 3, const int _m = 3) {
  int Nmax = nmax; // максимальное число итераций (не менее 1)
  std::vector<std::vector<double>> v; // сеточная функция ѵ (x, y)
  v.resize(_n + 1);
  for (size_t i = 0; i < _n + 1; i++) {
    v[i].resize(_m + 1);
    for (size_t j = 0; j < _m + 1; j++) {
      v[i][j] = 0;
    }
  }
  std::vector<std::vector<double>> v2; // сеточная функция ѵ (x, y)
  v2.resize(2 * _n + 1);
  for (size_t i = 0; i < 2 * _n + 1; i++) {
    v2[i].resize(2 * _m + 1);
    for (size_t j = 0; j < 2 * _m + 1; j++) {
      v2[i][j] = 0;
    }
  }
  std::vector<std::vector<double>> r_vec; // вектор невязки
  r_vec.resize(_n + 1);
  for (size_t i = 0; i < _n + 1; i++) {
    r_vec[i].resize(_m + 1);
    for (size_t j = 0; j < _m + 1; j++) {
      r_vec[i][j] = 0;
    }
  }
  std::vector<std::vector<double>> ar; // вектор невязки
  ar.resize(_n + 1);
  for (size_t i = 0; i < _n + 1; i++) {
    ar[i].resize(_m + 1);
    for (size_t j = 0; j < _m + 1; j++) {
      ar[i][j] = 0;
    }
  }
  double r1 = 0;
  int S = 0;         // счетчик итераций
  double eps = _eps; // заданная точность
  double eps_max = 0; // точность, достигнутая на текущей итерации
  double eps_cur = 0; // для подсчета точности на текущей итерации
  double a2, k2, h2; // ненулевые элементы матрицы (-А)
  const int n = _n, m = _m; // размерность сетки
  double r_max = 0;
  double a = 0;
  double b = 1;
  double c = 0;
  double d = 1; // границы области определения уравнения
  int i, j; // индексы
  double r = 0;
  double v_old; // старое значение преобразуемой компоненты вектора
  double v_new; // новое значение преобразуемой компоненты вектора и
  bool flag = false; // условие остановки
  double h = ((b - a) / n);
  double k = ((d - c) / m);
  double z = 0;
  double z_max = 0;
  double tau_s = 0;
  double num = 0; // tau=num/denominator
  double denominator = 1;
  h2 = -std::pow((n / (b - a)), 2);
  k2 = -std::pow((m / (d - c)), 2);
  a2 = -2 * (h2 + k2);
  for (int i = 0; i < n + 1; i++) {
    double x = a + i * h;
    v[i][0] = x - x * x;
    v[i][m] = x - x * x;
  }
  for (int j = 0; j < m + 1; j++) {
    double y = c + j * k;
    v[0][j] = std::sin(M_PI * y);
    v[n][j] = std::sin(M_PI * y);
  }
  do {
    z = 0;
    z_max = 0;
    eps_max = 0;
    r_max = -20;
    num = 0;
    denominator = 0;
    for (j = 1; j < m; j++) {
      for (i = n - 1; i > 0; i--) {
        r_vec[i][j] = f_main(a + i * h, c + j * k) -
                      (h2 * (v[i + 1][j] + v[i - 1][j]) +
                       k2 * (v[i][j + 1] + v[i][j - 1]) + v[i][j] * a2);
      }
    }
    for (j = 1; j < m; j++) {
      for (i = n - 1; i > 0; i--) {
        ar[i][j] =
            (h2 * (r_vec[i + 1][j] + r_vec[i - 1][j]) +
             k2 * (r_vec[i][j + 1] + r_vec[i][j - 1]) + r_vec[i][j] * a2);
        num += ar[i][j] * r_vec[i][j];
        denominator += ar[i][j] * ar[i][j];
      }
    }
    tau_s = num / denominator;
    for (j = 1; j < m; j++) {
      for (i = n - 1; i > 0; i--) {
        v_old = v[i][j];
        v_new = v_old + tau_s * r_vec[i][j];
        eps_cur = std::fabs(v_old - v_new);
        if (eps_cur > eps_max) {
          eps_max = eps_cur;
        }
        v[i][j] = v_new;
        r = fabs(f_main(a + i * h, c + j * k) -
                 (h2 * (v[i + 1][j] + v[i - 1][j]) +
                  k2 * (v[i][j + 1] + v[i][j - 1]) + v_old * a2));
        if (r > r_max)
          r_max = r;
      }
    }

    S = S + 1;
    if ((eps_max <= eps) or (S >= Nmax)) {
      flag = true;
    }
  } while (!flag);
  Write_main_param(n + 1, m + 1, S, eps_max, r_max);
  return v;
}
}

extern "C" {
int do_test(const int N, const int M, const double eps, const int steps) {
  int n = N;
  int m = M;
  std::vector<std::vector<double>> v(n, std::vector<double>(m, 0));
  MinimalResiduals_test(steps, eps, n, m);
  return 0;
}
}
extern "C" {
int do_main(const int N, const int M, const double eps, const int steps) {
  int n = N;
  int m = M;
  std::vector<std::vector<double>> v_1(n, std::vector<double>(m, 0));
  std::vector<std::vector<double>> v_2(n, std::vector<double>(2 * m, 0));
  v_1 = MinimalResiduals_main(steps, eps, n, m);
  v_2 = MinimalResiduals_main(steps, eps, 2 * n, 2 * m);
  n++;
  m++;
  double maxv_v2 = 0.0;
  double v_v2 = 0.0;
  int n2 = 0;
  int m2 = 0;
  for (int i = n - 1; i >= 0; i--) {
    for (int j = 0; j < m; j++) {
      v_v2 = std::fabs(v_1[i][j] - v_2[i * 2][j * 2]);
      if (v_v2 > maxv_v2) {
        maxv_v2 = v_v2;
        n2 = i;
        m2 = j;
      }
    }
  }
  Write_main_v(v_1, n, m, maxv_v2, n2, m2);
  return 0;
}
}
int main() {
  return 0;
}