// LLT разложение, матрица в профильном формате.
/*
* Структура программы:
* Ввод матрицы и вектора F
* Разложение матрицы на LLT составляющие
* Решение системы Ly = F  (прямой ход)
* Решение системы LTx = y (обратный ход)
*/

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

typedef float real;
typedef float real_acc;

using namespace std;

const string file1("matrix.txt"); // Имя файла, содержащего матрицу в профильном формате.
const string file2("vector.txt"); // Вектор - результат умножения матрицы на искомый вектор.
const string file3("result.txt"); // Искомый вектор - решение СЛАУ.
const string file4("subs.txt");   // Погрешность вычислений.

size_t N = 0; // Размерность матрицы
vector<size_t> ia;
vector<real> di;
vector<real> al;
vector<real> F;

vector<real> x;
vector<real> y; // Промежуточный вектор для прямого хода.

void input(string file1, string file2);
void LLT_dec();
void straight_stroke();
void reverse_stroke();
void cond_test_gen(size_t k, string file3, string file4);
void gilbert_test_gen(size_t k, string file3, string file4);
void Gauss_test_gen(size_t k, string file3, string file4);

void main()
{
   size_t k = 5;


   //input(file1, file2);
   //cond_test_gen(k, file3, file4);

   //gilbert_test_gen(k, file3, file4);

   Gauss_test_gen(k, file3, file4);





   N = 2;
   int SUMMA = 0;

   for (size_t k = 0; k < N; k++)
   {
      for (size_t i = 0; i < N + 1; i++)
         SUMMA++;

      for (size_t i = k + 1; i < N; i++)
      {
         SUMMA++;
         for (size_t j = 0; j < N + 1; j++)
            SUMMA += 2;
      }
   }

   cout << SUMMA;
   for (int k = N - 1; k > -1; k--)
   {
      for (int i = N; i > -1; i--)
         SUMMA++;
      for (int i = k - 1; i > -1; i--)
      {
         SUMMA++;
         for (int j = N; j > -1; j--)
            SUMMA += 2;
      }
   }
}

/*
* 1. Преобразование особых матрицы и вектора F
*    к нужному виду.
* 2. Решение СЛАУ.
* 3. Запись результата и погрешности относительно
*    точного решения.
*/
void cond_test_gen(size_t k, string file3, string file4)
{
   di[0] += 1.0 / pow(10, k);
   F[0]  += 1.0 / pow(10, k);

   LLT_dec();
   straight_stroke();
   reverse_stroke();

   // Вывод искомого вектора.
   size_t prec = 0;
   k < 7 ? prec = 7 : prec = 14;
   ofstream fout(file3);
   for (size_t i = 0; i < N; i++)
      fout << scientific << setprecision(prec) << x[i] << endl;
   fout.close();

   // Вывод погрешности относительно точного вектора.
   ofstream fout2(file4);
   vector<real> x_acc = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0 };
   for (size_t i = 0; i < 10; i++)
      fout2 << scientific << setprecision(2) << abs(x_acc[i] - x[i]) << endl;
   fout2.close();
}

/*
* 1. Создание матрицы Гильберта k-той размерности и
*    вектора F умножением на x(1, ..., N).
* 2. Решение СЛАУ.
* 3. Запись результата и погрешности относительно
*    точного решения.
*/
void gilbert_test_gen(size_t k, string file3, string file4)
{
   N = k;
   ia.resize(N + 1);
   ia[0] = 0;
   for (size_t i = 0; i < N; i++)
      ia[i + 1] = ia[i] + i;

   al.resize(ia[N]);
   di.resize(N);
   for (size_t i = 0; i < N; i++)
   {
      size_t i0 = ia[i];        // Индекс начального элемента i-той строки в массиве al.
      size_t i1 = ia[i + 1];    // Индекс начального элемента i + 1 стоки.
      size_t j = i - (i1 - i0); // Номер столбца первого элемента в строке.

      for (size_t ij = i0; ij < i1; ij++, j++)
         al[ij] = 1.0 / (i + j + 1);
      di[i] = 1.0 / (2 * i + 1);
   }

   F.resize(N);
   for (size_t i = 0; i < N; i++)
   {
      size_t ind = i + ia[i] - ia[i + 1];
      for (size_t j = ia[i]; j < ia[i + 1]; j++)
      {
         F[i] += al[j] * (ind + 1);
         F[ind] += al[j] * (i + 1);
         ind++;
      }
      F[i] += di[i] * (i + 1);
   }

   y.resize(N);
   x.resize(N);

   LLT_dec();
   straight_stroke();
   reverse_stroke();

   // Вывод искомого вектора.
   size_t prec = 0;
   k < 7 ? prec = 7 : prec = 14;
   ofstream fout(file3);
   fout << scientific << setprecision(14);
   for (size_t i = 0; i < N; i++)
      fout << x[i] << endl;
   fout.close();

   // Вывод погрешности относительно точного вектора.
   ofstream fout2(file4);
   vector<real> x_acc(N);
   for (size_t i = 0; i < N; i++)
      x_acc[i] = i + 1;
   fout2 << scientific << setprecision(2);
   for (size_t i = 0; i < N; i++)
      fout2 << abs(x_acc[i] - x[i]) << endl;
   fout2.close();
}

void input(string file1, string file2)
{
   // Переделать в сишные сканфс
   // для контроля количества вводимых знаков.

   ifstream fin(file1);
   ifstream fin2(file2);
   fin >> N;
   ia.resize(N + 1);
   for (size_t i = 0; i < N + 1; i++)
   {
      fin >> ia[i];
      ia[i]--;
   }

   di.resize(N);
   for (size_t i = 0; i < N; i++)
      fin >> di[i];

   al.resize(ia[N]);
   for (size_t i = 0; i < ia[N]; i++)
      fin >> al[i];
   fin.close();

   F.resize(N);
   for (size_t i = 0; i < N; i++)
      fin2 >> F[i];
   fin2.close();

   y.resize(N);
   x.resize(N);
}

// LLT - разложение.
// Заменяет матрицу в массивах di и al нижнетреугольной L
void LLT_dec()
{
   for (size_t i = 0; i < N; i++)
   {
      real_acc sumdi = 0.0;

      size_t i0 = ia[i];        // Индекс начального элемента i-той строки в массиве al.
      size_t i1 = ia[i + 1];    // Индекс начального элемента i + 1 стоки.
      size_t j = i - (i1 - i0); // Номер столбца первого элемента в строке.

      for (size_t ij = i0; ij < i1; ij++, j++)
      {
         real_acc suml = 0.0;

         size_t j0 = ia[j];       // Индекс начального элемента j-той строки в массиве al.
         size_t j1 = ia[j + 1];   // Индекс начального элемента j + 1 стоки.

         size_t ik = i0;
         size_t jk = j0;

         size_t ci = ij - i0;  // Количество элементов i-той строки.
         size_t cj = j1 - j0;  // Количество элементов j-той строки.

         // Сдвиг индексов начальных элементов строк
         // к первому существующему элементу наименьшей строки.
         ci > cj ? ik += ci - cj : jk += cj - ci;

         // Суммируем произведения элементов двух строк, стоящих до искомого.
         for (; ik < ij; ik++, jk++)
            suml += al[ik] * al[jk];
         al[ij] = (al[ij] - suml) / di[j];

         sumdi += al[ij] * al[ij];
      }
      di[i] = sqrt(di[i] - sumdi);
   }
}

// Прямой ход. (Ly = F)
// Заполнение вектора y.
void straight_stroke()
{
   for (size_t i = 0; i < N; i++)
   {
      size_t in = ia[i + 1];
      size_t in_prev = ia[i];

      size_t j = i - (in - in_prev);
      real_acc sum = 0.0;
      for (size_t k = in_prev; k < in; k++, j++)
         sum += al[k] * y[j];

      y[i] = (F[i] - sum) / di[i];
   }
}

// Обратный ход. (LTx = y)
// Заполнение искомого вектора.
void reverse_stroke()
{
   vector<real_acc> temp(N);
   fill(temp.begin(), temp.end(), 0.0);

   for (size_t i = N; i >= 1; i--)
   {
      size_t in = ia[i];
      size_t in_prev = ia[i - 1];

      x[i - 1] = (y[i - 1] - temp[i - 1]) / di[i - 1];
      size_t j = i - 2;
      for (size_t k = in; k > in_prev; k--, j--)
         temp[j] += al[k - 1] * x[i - 1];
   }
}

void Gauss_test_gen(size_t p, string file3, string file4)
{
   vector<vector<real_acc>> Matrix_Clone =
   {
      {  8.0,    0.0,  -3.0,   0.0,   0.0,   0.0,  -1.0,  -4.0,    0.0,    0.0,       -40.0 },
      {  0.0,    5.0,   0.0,  -1.0,   0.0,   0.0,   0.0,   0.0,   -2.0,   -2.0,       -32.0 },
      { -3.0,    0.0,   7.0,   0.0,  -1.0,   0.0,  -3.0,   0.0,    0.0,    0.0,       -8.0  },
      {  0.0,   -1.0,   0.0,   8.0,   0.0,  -2.0,   0.0,   0.0,   -1.0,   -4.0,       -31.0 },
      {  0.0,    0.0,  -1.0,   0.0,   4.0,   0.0,   0.0,  -3.0,    0.0,    0.0,       -7.0  },
      {  0.0,    0.0,   0.0,  -2.0,   0.0,   6.0,  -3.0,   0.0,   -1.0,    0.0,       -2.0  },
      { -1.0,    0.0,  -3.0,   0.0,   0.0,  -3.0,   9.0,   0.0,    0.0,   -2.0,        15.0 },
      { -4.0,    0.0,   0.0,   0.0,  -3.0,   0.0,   0.0,   11.0,  -4.0,    0.0,        33.0 },
      {  0.0,   -2.0,   0.0,  -1.0,   0.0,  -1.0,   0.0,  -4.0,    11.0,  -3.0,        23.0 },
      {  0.0,   -2.0,   0.0,  -4.0,   0.0,   0.0,  -2.0,   0.0,   -3.0,    11.0,       49.0 }
   };

   N = Matrix_Clone.size();
   Matrix_Clone[0][0] += 1.0 / pow(10, p);
   Matrix_Clone[0][N] += 1.0 / pow(10, p);

   vector<vector <real>> Matrix(N, vector <real>(N + 1));
   for (size_t i = 0; i < N; i++)
      for (size_t j = 0; j < N + 1; j++)
         Matrix[i][j] = Matrix_Clone[i][j];

   // Прямой ход (Зануление нижнего левого угла)
   for (size_t k = 0; k < N; k++)      // k-номер строки.
   {
      for (size_t i = 0; i < N + 1; i++)      // i-номер столбца.
         Matrix_Clone[k][i] /= Matrix[k][k];  // Деление k-строки на первый член !=0 для преобразования его в единицу.

      for (size_t i = k + 1; i < N; i++)                         // i-номер следующей строки после k.
      {
         real K = Matrix_Clone[i][k] / Matrix_Clone[k][k];       // Коэффициент.
         for (size_t j = 0; j < N + 1; j++)                      // j-номер столбца следующей строки после k.
            Matrix_Clone[i][j] -= Matrix_Clone[k][j] * K;        // Зануление элементов матрицы ниже первого члена, преобразованного в единицу
      }

      for (size_t i = 0; i < N; i++)                             // Обновление, внесение изменений в начальную матрицу
         for (size_t j = 0; j < N + 1; j++)
            Matrix[i][j] = Matrix_Clone[i][j];
   }

   // Обратный ход (Зануление верхнего правого угла)
   for (int k = N - 1; k > -1; k--) // k-номер строки
   {
      for (int i = N; i > -1; i--)  // i-номер столбца
         Matrix_Clone[k][i] /= Matrix[k][k];
      for (int i = k - 1; i > -1; i--) //i-номер следующей строки после k
      {
         real K = Matrix_Clone[i][k] / Matrix_Clone[k][k];
         for (int j = N; j > -1; j--) //j-номер столбца следующей строки после k
            Matrix_Clone[i][j] -= Matrix_Clone[k][j] * K;
      }
   }

   // Отделяем от общей матрицы ответы
   vector<real> x(N);
   for (size_t i = 0; i < N; i++)
      x[i] = Matrix_Clone[i][N];

   // Вывод искомого вектора.
   size_t prec = 0;
   p < 7 ? prec = 7 : prec = 14;
   ofstream fout(file3);
   for (size_t i = 0; i < N; i++)
      fout << scientific << setprecision(prec) << x[i] << endl;
   fout.close();

   // Вывод погрешности относительно точного вектора.
   ofstream fout2(file4);
   vector<real> x_acc = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0 };
   for (size_t i = 0; i < 10; i++)
      fout2 << scientific << setprecision(2) << abs(x_acc[i] - x[i]) << endl;
   fout2.close();
}