ChM_lab1

#include <cstdio>
#include <cstring>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <vector>
#include <string>
#include <sstream>
#include <cmath>
#define EPS 1E-14

double dabs(double _) { return _ > 0 ? _ : -_; }

// Матрица в профильном формате
template<typename T>
struct SkylineStorageMatrix
{
    std::vector<int> ia; // индексы начала каждого столбца
    std::vector<T> al;   // элементы ниже глав. диагонали
    std::vector<T> au;   // элементы выше глав. диагонали
    std::vector<T> di;   // элементы глав. диагонали
};

// Считать матрицу в профильном формате
template<typename T>
void readSkyline(std::istream& from, SkylineStorageMatrix<T>& to)
{
    int n;
    from >> n;                        // размер матрицы
    to.ia.resize(n + 1);              // резерв места для индексов (n+1)
    for (int i = 0; i < n + 1; i++)
    {
        from >> to.ia[i];         // индексы начала каждого столбца
    }

    int k = to.ia[n];             // количество ненулевых элементов
    to.al.resize(k);              // резерв места для элементов нижнего треугольника
    for (int i = 0; i < k; i++)
    {
        from >> to.al[i];         // записываем элементы ниже глав. диагонали
    }
    to.au.resize(k);              // резерв места для элементов выше глав. диагонали
    for (int i = 0; i < k; i++)
    {
        from >> to.au[i];         // записываем элементы выше глав. диагонали
    }

    to.di.resize(n);
    for (int i = 0; i < n; i++)   // резерв места для элементов глав. диагонали
    {
        from >> to.di[i];         // записываем элементы глав. диагонали
    }
}

// Сконвертировать матрицу в профильном формате в квадратную
template<typename T>
void skylineToSquare(SkylineStorageMatrix<T>& a, std::vector<std::vector<T>>& s)
{
    // вычисление кол-ва ненулевых элементов в строке i
    #define elementsInLine(i) (a.ia[i + 1] - a.ia[i])

    // изменяем размер матрицы в соотв. с размером главной диагонали
    s.resize(a.di.size());

    // заполнение матрицы нулями
    for (unsigned int i = 0; i < a.di.size(); i++)
    {
        s[i].resize(a.di.size());   // устанавливаем размер каждой строки
        for (unsigned int j = 0; j < a.di.size(); j++)
        {
            s[i][j] = (T)0;         // каждый элемент приравниваем к 0
        }
    }

    // заполнение главной диагонали матрицы значениями из вектора
    for (unsigned int i = 0; i < a.di.size(); i++)
    {
        s[i][i] = a.di[i];
    }

    // заполнение элементов ниже и выше глав. диагонали
    for (unsigned int i = 1; i < a.di.size(); i++)
    {
        // считаем смещение для заполнения элементов ниже глав. диагонали
        int offset = i - elementsInLine(i);
        for (int j = 0; j < elementsInLine(i); j++)
        {
            s[i][offset + j] = a.al[a.ia[i] + j];  // заполнение элементов ниже глав. диагонали
                                                   // это все классно, заполняем симметричные
            s[offset + j][i] = a.au[a.ia[i] + j];  // заполнение элементов выше глав. диагонали
        }
    }
    // больше не считаем ненулевые элементы
    #undef elementsInLine
}

 // Считать вектор
template<typename T>
void readVector(std::istream& from, std::vector<T>& to)
{
    int n;         // берем размерность вектора
    from >> n;
    to.resize(n);  // изменяем размер вектора to
    for (int i = 0; i < n; i++)
    {
        from >> to[i]; // передаем параметры с ввода в ячейки массива
    }
}

// Распечатать квадратную матрицу
template<typename T>
void printMatrix(const std::vector<std::vector<T>> a)
{
    for (unsigned int i = 0; i < a.size(); i++)  // обычный вывод таблицы
    {
        for (unsigned int j = 0; j < a.size(); j++)
        {
            std::cout << std::setw(20) << a[i][j];
        }
        std::cout << std::endl;
    }
}

// Распечатать матрицу в профильном формате
template<typename T>
void printSkyline(SkylineStorageMatrix<T>& a)
{
    std::vector<std::vector<T>> s;
    skylineToSquare(a, s); // берем 4 вектора, выводим их в разные строки.
    printMatrix(s);        // выводим как матрицу
}


// LU(sq)-разложение матрицы в профильном формате
template<typename T, typename T2=T>                // создание двух типов для быстрого переключения
bool LusqSkyline(SkylineStorageMatrix<T>& a)
{
    // Присвоить dst значение корня из val.
    // Если val <= 0, вернуть false из текущей подпрограммы.
    // проверка на отрицательный корень
    #define _sqrt(dst, val) { double underSqrt = (double)(val); if (underSqrt <= 0) return false; dst = (T)sqrt(underSqrt); }

    // Количество элементов в i'ой строке/столбце (ненулевых элементов)
    #define elementsInLine(i) (a.ia[i + 1] - a.ia[i])

    // Первая строка матрицы
    _sqrt(a.di[0], a.di[0]); // d_1 = sqrt(a_11);
                             // вычисляется квадратный корень из первого элемента диагонали
                             // Если этот элемент меньше или равен нулю, функция завершится с возвратом false.
    for (unsigned int i = 1; i < a.di.size(); i++)
    {   // Нужны ненулевые элементы первой строки
        // i'ый элемент первой строки ненулевой <=> i'ый столбец содержит i элементов
        if (elementsInLine(i) == i)
        {
            a.au[a.ia[i]] = a.au[a.ia[i]] / a.di[0]; // U_1i = a_1i/q1
            // a.au[a.ia[i]] - что это.
        }
    }

    // Обход по строкам матрицы
    for (unsigned int i = 1; i < a.di.size(); i++)
    {
        if (a.al.size())
        {
            // l_ij = (a_ij - \sum_{k=1}{j-1} l_ik * u_kj) / d_j
            // Для вычисления j'го элемента i'ой строки, необходимо найти
            // сумму произведений первых j-1 элементов i'ой строки на
            // соответствующие первые j-1 элементов j'го столбца.
            // Затем вычесть её из a_ij и полученную разность поделить на a_jj.
            // "Легко заметить, что подобное преобразование сохраняет портрет матрицы"

            unsigned int alStartIndex = a.ia[i];                   // индекс первого ненулевого элемента i'ой строки в массиве al
            unsigned int alEndIndex = a.ia[i] + elementsInLine(i); // индекс последнего ненулевого элемента (????)
            unsigned int j = i - elementsInLine(i);                // индекс текущего столбца

            // До первого ненулевого элемента не было ненулевых элементов (??????)
            a.al[alStartIndex] = a.al[alStartIndex] / a.di[j];     // нормализация первого элемента
            alStartIndex++;                                        // переход к сл. элементу строки
            j++;                                                   // переход к сл. столбцу

            // Основной цикл
            for (unsigned int alIndex = alStartIndex; alIndex < alEndIndex; alIndex++, j++)
            {
                T sum = (T)0;                                      // Инициализация суммы
                unsigned int elementsInLineI = elementsInLine(i);  // Количество ненулевых элементов в строке i
                unsigned int elementsInColJ = elementsInLine(j);   // Количество ненулевых элементов в столбце j

                // Для каждого ненулевого элемента из [alIndex...alEndIndex]
                // вычисляем сумму произведений элементов l_ik и u_kj
                // l_ik — элемент нижней треугольной матрицы, соответствующий текущему элементу строки.
                // u_kj — элемент верхней треугольной матрицы, соответствующий текущему элементу столбца.
                // Если индекс k меньше количества ненулевых элементов
                // в текущем столбце или строке, то получаем соответствующие значения,
                // иначе присваиваем 0.

                // Цикл для вычисления суммы произведений
                for (unsigned int k = 0; k < j; k++)
                {
                    T l_ik = (k < i - elementsInLineI) ? (T)0 : a.al[a.ia[i] + (k - (i - elementsInLineI))];
                    T u_kj = (k < j - elementsInColJ) ? (T)0 : a.au[a.ia[j] + (k - (j - elementsInColJ))];
                    sum += l_ik * u_kj;
                }
                // Обновляем элемент alIndex
                a.al[alIndex] = (a.al[alIndex] - sum) / a.di[j];
            }
        }

        // d_i = sqrt(d_i - \sum_{k=1}{i-1} l_ik * u_ki) // отладка
        // Повторение основного цикла для T2

        {
            T2 sum = (T2)0;
            unsigned int elementsInLineI = elementsInLine(i); // Получаем количество ненулевых элементов в строке i
            for (unsigned int k = 0; k < i; k++)
            {
                T l_ik = (k < i - elementsInLineI) ? (T)0 : a.al[a.ia[i] + (k - (i - elementsInLineI))];
                T l_ki = (k < i - elementsInLineI) ? (T)0 : a.au[a.ia[i] + (k - (i - elementsInLineI))];
                sum += (T2)l_ik * (T2)l_ki;  // Приведение к типу T2 для суммирования
            }
            // printf("%d %lf\n", i, sum); // отладка
            _sqrt(a.di[i], a.di[i] - sum); // обновляем диагональный элемент
        }

        if (a.au.size())  // проверка на наличие элементов в верхней матрице
        {
            // u_ij = (a_ij - \sum_{k=1}{i-1} l_ik * u_kj) / d_j
            // То же самое, что l_ij, только писать нужно в al (не последовательно) и сумма по k до i,
            // и делить на di[i]

            for (unsigned int j = i + 1; j < a.di.size(); j++)
            {
                if (i >= j - elementsInLine(j)) // Проверка, что текущий индекс i не выходит за пределы j
                                                // т.к матрица квадратная
                {
                    T2 sum = (T2)0;
                    unsigned int elementsInLineI = elementsInLine(i); // количество ненулевых элементов в строке
                    unsigned int elementsInColJ = elementsInLine(j); // количество ненулевых элементов в столбце
                    for (unsigned int k = 0; k < i; k++)
                    {
                        T l_ik = (k < i - elementsInLineI) ? (T)0 : a.al[a.ia[i] + (k - (i - elementsInLineI))];
                        T u_kj = (k < j - elementsInColJ) ? (T)0 : a.au[a.ia[j] + (k - (j - elementsInColJ))];
                        sum += (T2)l_ik * (T2)u_kj;
                    }

                    unsigned int auIndex = a.ia[j] + (i - (j - elementsInLine(j))); // Вычислим индекс для записи в верхнюю треугольную матрицу
                    a.au[auIndex] = (a.au[auIndex] - sum) / a.di[i]; // Обновляем значение в верхней треугольной матрице
                }
            }
        }
    }

    #undef elementsInLine
    #undef _sqrt
    // printSkyline(a); // отладка
    return true;
}

// Решить уравнение Ly=b
template<typename T, typename T2=T>
bool findY(std::vector<T>& y, SkylineStorageMatrix<T>& a, std::vector<T>& b)
{
    // dst = частное a и b.
    // Если b = 0, вернуть false из текущей подпрограммы.

    #define _div(dst, a, b) { if (((b) > (T)0 ? (b) : -(b)) < (T)EPS) return false; dst = (a) / (b); }

    // Количество элементов в i'ой строке/столбце
    #define elementsInLine(i) (a.ia[i + 1] - a.ia[i])

    for (unsigned int i = 0; i < a.di.size(); i++)
    {
        T2 sum = (T2)0;
        unsigned int elementsInLineI = elementsInLine(i); // кол-во ненулевых элементов
        // цикл по ним
        for (unsigned int alIndex = a.ia[i], bIndex = i - elementsInLineI, counter = 0; counter < elementsInLineI; alIndex++, bIndex++, counter++)
        {
            sum += (T2)a.al[alIndex] * (T2)y[bIndex]; // суммирование
        }

        _div(y[i], b[i] - sum, a.di[i]);         // Делим (b[i] - sum) на a.di[i] и сохраняем результат в y[i]
    }

    #undef _div
    #undef elementsInLine
    return true;
}

// Решить уравнение Ux=y
template<typename T, typename T2=T>
bool findX(std::vector<T>& x, SkylineStorageMatrix<T>& a, std::vector<T>& y)
{
    #define _div(dst, a, b) { if (((b) > (T)0 ? (b) : -(b)) < (T)EPS) return false; dst = (a) / (b); }

    // Количество элементов в i'ой строке/столбце
    #define elementsInLine(i) (a.ia[i + 1] - a.ia[i])

    for (unsigned int i = a.di.size() - 1; true; i--)
    {
        T2 sum = (T2)0;
        for (unsigned int j = i + 1; j < a.di.size(); j++)
        {
            unsigned int elementsInColJ = elementsInLine(j);
            // Проверяем, есть ли ненулевые элементы в колонне j, которые влияют на текущую строку i
            sum += (i < j - elementsInColJ) ? (T2)0 : (T2)a.au[a.ia[j] + (i - (j - elementsInColJ))] * x[j];
        }
        // Делим (y[i] - sum) на a.di[i] и сохраняем результат в x[i]
        _div(x[i], y[i] - sum, a.di[i]);
        if (i == 0) break; // Прерываем цикл, если достигли первой строки
    }

    #undef _div
    #undef elementsInLine
    return true;
}

// Решить уравнение Ax=B
template<typename T, typename T2=T>
bool solve(std::vector<T>& x, SkylineStorageMatrix<T>& a, std::vector<T>& b)
{
    std::vector<T> y(b.size()); // вектор "y" с размером "b"

    // Проверка на возможность LU-разложения
    if (!LusqSkyline<T,T2>(a))
    {
        std::cerr << "Matrix is not LU(sq)-decomposeable" << std::endl;
        return false;
    }
    // Решение промежуточной системы y = L * b
    // методом проверки "L * y = b" и "U * x = y"
    if (!findY<T,T2>(y, a, b) || !findX<T,T2>(x, a, y))
    {
        std::cerr << "System is inconsistent" << std::endl;
        return false;
    }

    return true;
}

// Записать матрицу в профильном формате
template<typename T>
void writeSkyline(SkylineStorageMatrix<T>& from, std::ostream& to)
{
    to << from.di.size() << std::endl; // запись размера диагональных элементов
    for (unsigned int i = 0; i < from.ia.size(); i++)
    {
        to << from.ia[i] << " ";
    }
    to << std::endl;

    for (unsigned int i = 0; i < from.al.size(); i++) // запись массива al
    {
        to << std::fixed << std::setprecision(14) << from.al[i] << " ";
    }
    to << std::endl;
    for (unsigned int i = 0; i < from.au.size(); i++) // au
    {
        to << std::fixed << std::setprecision(14) << from.au[i] << " ";
    }
    to << std::endl;

    for (unsigned int i = 0; i < from.di.size(); i++) // di
    {
        to << std::fixed << std::setprecision(14) << from.di[i] << " ";
    }
    to << std::endl;
}

// Записать вектор
template<typename T>
void writeVector(std::vector<T>& from, std::ostream& to)
{
    to << from.size() << std::endl;  // запись размера массива

    for (unsigned int i = 0; i < from.size(); i++)
    {
        to << std::fixed << std::setprecision(14) << from[i] << std::endl;
        // каждый элемент выводится в фиксированном формате с точностью 14 знаков
    }
    to << std::endl;
}

// Провести серию тестов на число обусловленности
void conditionNumberTestSeries(int tests)
{
    // создание таблицы для вывода
    std::cout << "\\newline{\\tiny\\tt\\begin{tabular}{| l | p{3cm} | p{1cm} | p{3cm} | p{1cm} | p{3cm} | p{1cm} |}\n\\hline\n";
    std::cout << "k & $x^k$ (float) & $x^* - x^k$ (float) & $x^k$ (double) & $x^* - x^k$ (double) & $x^k$ (mixed) & $x^* - x^k$ (mixed) \\\\ \\hline\n";

    for (int test = 0; test < tests; test++)
    {
        // чтение матрица и вектора float
        std::fstream f_ifile("test6.txt");
        struct SkylineStorageMatrix<float> f_a;
        readSkyline(f_ifile, f_a);              // матрица
        int n = f_a.di.size();
        std::vector<float> f_b;
        readVector(f_ifile, f_b);               // вектор
        f_a.di[0] += (float)pow(10.0, -test);   // проверка на устойчивость
        f_b[0] += (float)(pow(10.0, -test));
        // в первый элемент диагонального вектора матрицы и вектора b добавляется небольшое значение
        // зависимое от номера теста
        std::vector<float> f_x(n);
        bool ok1 = solve(f_x, f_a, f_b);

        // double
        std::fstream d_ifile("test6.txt");
        struct SkylineStorageMatrix<double> d_a;
        readSkyline(d_ifile, d_a);
        std::vector<double> d_b;
        readVector(d_ifile, d_b);
        d_a.di[0] += (double)pow(10.0, -test);
        d_b[0] += (double)(pow(10.0, -test));
        std::vector<double> d_x(n);
        bool ok2 = solve(d_x, d_a, d_b);

        // mixed
        std::fstream m_ifile("test6.txt");
        struct SkylineStorageMatrix<float> m_a;
        readSkyline(m_ifile, m_a);
        std::vector<float> m_b;
        readVector(m_ifile, m_b);
        m_a.di[0] += (float)pow(10.0, -test);
        m_b[0] += (float)(pow(10.0, -test));
        std::vector<float> m_x(n);
        bool ok3 = solve<float,double>(m_x, m_a, m_b);

        // std::scientific, std::fixed - добавляют точность при форматировании
        //

        if (ok1 || ok2)
        {
            std::cout << test << " & ";
            for (int i = 0; i < n; i++)
            {
                std::cout<< std::fixed << std::setprecision(7) << f_x[i] << "\\newline";
            }
            std::cout << " & ";
            for (int i = 0; i < n; i++)
            {
                std::cout << std::scientific << std::setprecision(2) << fabs((i + 1) - f_x[i]) << "\\newline";
            }
            std::cout << " & ";
            for (int i = 0; i < n; i++)
            {
                std::cout << std::fixed << std::setprecision(14) << d_x[i] << "\\newline";
            }
            std::cout << " & ";
            for (int i = 0; i < n; i++)
            {
                std::cout << std::scientific << std::setprecision(2) << dabs((i + 1) - d_x[i]) << "\\newline";
            }
            std::cout << " & ";
            for (int i = 0; i < n; i++)
            {
                std::cout << std::fixed << std::setprecision(7) << m_x[i] << "\\newline";
            }
            std::cout << " & ";
            for (int i = 0; i < n; i++)
            {
                std::cout << std::scientific << std::setprecision(2) << fabs((i + 1) - m_x[i]) << "\\newline";
            }
            std::cout << " \\\\ \\hline \n";
        }

        if (test != tests - 1)
        {
            std::cout << "\\end{tabular}\\newline\\newline\\newline\n";
            std::cout << "\\begin{tabular}{| l | p{3cm} | p{1cm} | p{3cm} | p{1cm} | p{3cm} | p{1cm} |}\n\\hline\n";
            std::cout << "k & $x^k$ (float) & $x^* - x^k$ (float) & $x^k$ (double) & $x^* - x^k$ (double) & $x^k$ (mixed) & $x^* - x^k$ (mixed) \\\\ \\hline\n";
        }
    }

    std::cout << "\\end{tabular}}\n";
}


// Решить систему Ax=B с плотной матрицей методом Гаусса с выбором ведущего элемента
template<typename T>
bool solveGauss(std::vector<std::vector<T>> A, std::vector<T>& x, std::vector<T> B)
{
    unsigned int n = A.size();
    std::vector<std::vector<T>> t = A;

    for (unsigned int k = 0; k < n - 1; k++)
    {
        // Постолбцовый выбор главного элемента
        int m = k;
        for (int i = k + 1; i < n; i++)
        {
            if (abs(A[i][k]) > abs(A[m][k]))
            {
                m = i;
            }
        }
        if (abs(A[m][k]) < 1E-5)
        {
            return false; // Система не имеет уникального решения
        }
        for (int j = k; j < n; j++)
        {
            T tmp = A[k][j];
            A[k][j] = A[m][j];
            A[m][j] = tmp;
        }
        T tmp = B[k];
        B[k] = B[m];
        B[m] = tmp;
        // Конец

        for (int i = k + 1; i < n; i++)
        {
            t[i][k] = A[i][k] / A[k][k];
            B[i] = B[i] - t[i][k] * B[k];

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

    x[n - 1] = B[n - 1] / A[n - 1][n - 1];
    for (int k = n - 2; k >= 0; k--)
    {
        T sum = 0.0;
        for (int j = k + 1; j < n; j++)
        {
            sum += A[k][j] * x[j];
        }

        x[k] = (B[k] - sum) / A[k][k];
    }
}

// Сравнить Гаусса и разложение тестами на числом обусловленности
void compareGaussLusq(int tests)
{
    std::cout << "\\newline{\\tiny\\tt\\begin{tabular}{| l | p{3cm} | p{1cm} | p{3cm} | p{1cm} |}\n\\hline\n";
    std::cout << "k & $x^k$ (LU(sq)) & $x^* - x^k$ (LU(sq)) & $x^k$ (Gauss) & $x^* - x^k$ (Gauss) \\\\ \\hline\n";

    for (int test = 0; test < tests; test++)
    {
        std::fstream d_ifile("test6.txt");
        struct SkylineStorageMatrix<double> d_a;
        readSkyline(d_ifile, d_a);
        int n = d_a.di.size();
        std::vector<double> d_b;
        readVector(d_ifile, d_b);
        d_a.di[0] += (double)pow(10.0, -test);
        d_b[0] += (double)(pow(10.0, -test));
        std::vector<double> d_x(n);
        bool ok1 = solve(d_x, d_a, d_b);

        std::fstream g_ifile("test6.txt");
        struct SkylineStorageMatrix<double> g_a;
        readSkyline(g_ifile, g_a);
        std::vector<double> g_b;
        readVector(g_ifile, g_b);
        g_a.di[0] += (double)pow(10.0, -test);
        g_b[0] += (double)(pow(10.0, -test));
        std::vector<std::vector<double>> A;
        skylineToSquare(g_a, A);
        std::vector<double> g_x(n);
        bool ok2 = solveGauss(A, g_x, g_b);

        if (ok1 || ok2)
        {
            std::cout << test << " & ";
            for (int i = 0; i < n; i++)
            {
                std::cout<< std::fixed << std::setprecision(14) << d_x[i] << "\\newline";
            }
            std::cout << " & ";
            for (int i = 0; i < n; i++)
            {
                std::cout << std::scientific << std::setprecision(2) << dabs((i + 1) - d_x[i]) << "\\newline";
            }
            std::cout << " & ";
            for (int i = 0; i < n; i++)
            {
                std::cout << std::fixed << std::setprecision(14) << g_x[i] << "\\newline";
            }
            std::cout << " & ";
            for (int i = 0; i < n; i++)
            {
                std::cout << std::scientific << std::setprecision(2) << dabs((i + 1) - g_x[i]) << "\\newline";
            }
            std::cout << " \\\\ \\hline \n";
        }

        if (test != tests - 1)
        {
            std::cout << "\\end{tabular}\\newline\\newline\\newline\n";
            std::cout << "\\begin{tabular}{| l | p{3cm} | p{1cm} | p{3cm} | p{1cm} |}\n\\hline\n";
            std::cout << "k & $x^k$ (LU(sq)) & $x^* - x^k$ (LU(sq)) & $x^k$ (Gauss) & $x^* - x^k$ (Gauss) \\\\ \\hline\n";
        }
    }

    std::cout << "\\end{tabular}}\n";
}

// Сгенерировать матрицу Гильберта
template<typename T>
void generateHilbert(SkylineStorageMatrix<T>& dst, std::vector<T>& dstr, int n)
{
    // Создаем квадратную матрицу для удобства заполнения значениями
    std::vector<std::vector<T>> square(n);
    for (int i = 0; i < n; i++)
    {
        square[i].resize(n);    // Инициализируем каждую строку квадратной матрицы
    }

    // Инициализация массива индексов начала строк для профильной матрицы
    dst.ia.resize(n + 1);
    dst.ia[0] = 0;
    dst.ia[1] = 0;
    for (int i = 2; i < n + 1; i++)
    {
        dst.ia[i] = dst.ia[i - 1] + (i - 1);   // Количество ненулевых элементов в строке i равно i - 1
    }

    int ii, jj;     // текущая строка, столбец
    T sumij = (T)0;

    // Инициализируем массив для хранения ненулевых элементов в профильном формате
    dst.al.resize(dst.ia[n]);
    ii = 1, jj = 0; // 2 строка, 1 столбец
    for (int i = 0; i < dst.ia[n]; i++)
    {
        // Заполняем элемент матрицы Гильберта по формуле H_{ij} = 1 / (i + j + 1)
        dst.al[i] = (T)((T)1 / (T)(ii + jj + 1));
        // Сохраняем элемент в квадратную матрицу для удобства
        sumij += dst.al[i];

        square[ii][jj] = dst.al[i];         // Сохраняем элемент в квадратную матрицу для удобства
        if (++jj == ii)                     // Переход к следующему элементу в строке
        {
            // Если достигли конца текущей строки, переходим на следующую
            // Увеличиваем номер строки и сбрасываем номер столбца
            ii++, jj = 0;
        }
    }

    // Инициализация массива для хранения
    // ненулевых элементов в профильном формате
    dst.au.resize(dst.ia[n]); // резервируем память
    ii = 1, jj = 0;           // 1 столбец 2 строка
    for (int i = 0; i < dst.ia[n]; i++)
    {
        dst.au[i] = (T)((T)1 / (T)(ii + jj + 1));
        sumij += dst.au[i]; /// куда ведет переменная(???)

        square[jj][ii] = dst.al[i];
        if (++jj == ii)
        {
            ii++, jj = 0;
        }
    }

    dst.di.resize(n);
    for (int i = 0; i < n; i++)
    {
        dst.di[i] = (T)((T)1 / (T)(i + i + 1));;
        square[i][i] = dst.di[i]; // заполняем диагональные по формуле Гильберта
    }

    // Вектор
    dstr.resize(n);
    for (int i = 0; i < n; i++)
    {
        dstr[i] = (T)0;
        for (int j = 0; j < n; j++)
        {
            // Вычисляем элемент dstr[i] как сумму произведений элементов матрицы square и индексов
            dstr[i] += square[i][j] * (T)(j + 1);
        }
    }

    // printMatrix(square); // отладка
    std::ostringstream os; os << n;
    std::fstream out("hilbert" + std::string(os.str()) + ".txt", std::ios::out);
    writeSkyline(dst, out);  // Записываем профильную матрицу в файл
}

// Провести серию тестов с матрицами Гильберта размерности от n1 до n2
void hilbertTestSeries(int n1, int n2)
{
    std::cout << "\\newline{\\tiny\\tt\\begin{tabular}{| l | p{3cm} | p{3cm} | p{3cm} | p{3.1cm} |}\n\\hline\n";
    std::cout << "k & $x^k$ (float) & $x^* - x^k$ (float) & $x^k$ (double) & $x^* - x^k$ (double) \\\\ \\hline\n";

    for (int n = n1; n <= n2; n++)
    {
        struct SkylineStorageMatrix<float> f_a;
        std::vector<float> f_b;
        generateHilbert(f_a, f_b, n);
        std::vector<float> f_x(n);
        bool ok1 = solve(f_x, f_a, f_b);

        struct SkylineStorageMatrix<double> d_a;
        std::vector<double> d_b;
        generateHilbert(d_a, d_b, n);
        std::vector<double> d_x(n);
        bool ok2 = solve(d_x, d_a, d_b);

        std::cout << n << " & ";
        for (int i = 0; i < n; i++)
        {
            std::cout << std::fixed << std::setprecision(7) << (ok1 ? f_x[i] : 0.0) << "\\newline";
        }
        std::cout << " & ";
        for (int i = 0; i < n; i++)
        {
            std::cout << std::scientific << std::setprecision(2) << fabs(ok1 ? ((i + 1) - f_x[i]) : 0.0) << "\\newline";
        }
        std::cout << " & ";
        for (int i = 0; i < n; i++)
        {
            std::cout << std::fixed << std::setprecision(14) << (ok2 ? d_x[i] : 0.0) << "\\newline";
        }
        std::cout << " & ";
        for (int i = 0; i < n; i++)
        {
            std::cout << std::scientific << std::setprecision(2) << dabs(ok2 ? ((i + 1) - d_x[i]) : 0.0) << "\\newline";
        }
        std::cout << " \\\\ \\hline \n";

        if (n != n2)
        {
            std::cout << "\\end{tabular}\\newline\\newline\\newline\n";
            std::cout << "\\begin{tabular}{| l | p{3cm} | p{3cm} | p{3cm} | p{3.1cm} |}\n\\hline\n";
            std::cout << "k & $x^k$ (float) & $x^* - x^k$ (float) & $x^k$ (double) & $x^* - x^k$ (double) \\\\ \\hline\n";
        }
    }

    std::cout << "\\end{tabular}}\n";
}

int main(int argc, char** argv)
{
     // Проверяем, что передано
     // 4 аргумента и первый аргумент - "solve-double"
    if (argc == 4 && std::string(argv[1]) == "solve-double")
    {
        std::fstream ifile(argv[2], std::ios::in); // Открываем входной файл для чтения
        struct SkylineStorageMatrix<double> a;// Создаем объект для хранения матрицы в профильном формате
        readSkyline(ifile, a);  // Читаем профильную матрицу из файла
        std::vector<double> b; // Вектор для хранения правой части системы
        readVector(ifile, b);  // Читаем вектор из файла

        std::fstream ofile(argv[3], std::ios::out); // Открываем выходной файл для записи
        std::vector<double> x(b.size()); // Вектор для хранения решения
        solve(x, a, b);  // Решаем систему уравнений
        writeVector(x, ofile); // Записываем решение в выходной файл
    }
    else if (argc == 4 && std::string(argv[1]) == "solve-float")
    {
        std::fstream ifile(argv[2], std::ios::in);
        struct SkylineStorageMatrix<float> a;
        readSkyline(ifile, a);
        std::vector<float> b;
        readVector(ifile, b);

        std::fstream ofile(argv[3], std::ios::out);
        std::vector<float> x(b.size());
        solve(x, a, b);
        writeVector(x, ofile);
    }
    else if (argc == 4 && std::string(argv[1]) == "solve-double-square")
    {
        std::fstream ifile(argv[2], std::ios::in);
        struct SkylineStorageMatrix<double> a;
        readSkyline(ifile, a);
        std::vector<double> b;
        readVector(ifile, b);

        std::vector<std::vector<double>> A;
        skylineToSquare(a, A);

        std::fstream ofile(argv[3], std::ios::out);
        std::vector<double> x(b.size());
        solveGauss(A, x, b); // Решаем систему уравнений методом Гаусса
        writeVector(x, ofile);
    }
    else if (argc == 4 && std::string(argv[1]) == "hilbert")
    {
        std::istringstream os1(argv[2]); int from; os1 >> from;
        std::istringstream os2(argv[3]); int to; os2 >> to;
        hilbertTestSeries(from, to); // Запускаем тестирование матриц Гильберта
    }
    else if (argc == 3 && std::string(argv[1]) == "condition")
    {
        std::istringstream os(argv[2]); int n; os >> n;
        conditionNumberTestSeries(18); // Запускаем тестирование числа обусловленности
    }
    else if (argc == 3 && std::string(argv[1]) == "competition")
    {
        std::istringstream os(argv[2]); int n; os >> n;
        compareGaussLusq(18); // Сравниваем методы Гаусса и LU-разложения
    }
    else
    {
        // Выводим сообщение о правильном использовании программы
        std::cerr << "Usage: " << argv[0] << " solve-double test.txt output.txt" << std::endl
                  << "       " << argv[0] << " solve-float test.txt output.txt" << std::endl
                  << "       " << argv[0] << " solve-double-square test.txt output.txt" << std::endl
                  << "       " << argv[0] << " solve-float-square test.txt output.txt" << std::endl
                  << "       " << argv[0] << " hilbert from to" << std::endl
                  << "       " << argv[0] << " condition n" << std::endl
                  << "       " << argv[0] << " competition n" << std::endl;
    }

    return 0;
}