tridiag.cpp

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

using namespace std;

ifstream fin ("in.txt");
ofstream fout ("out_tridiag.txt");

const double eps = 10e-9;

void print(const vector<vector<double>> & v) {
    for (auto & row : v) {
        for (auto & it : row)
            fout << setw(12) << (fabs(it) > eps ? it : 0);
        fout << endl;
    }
}

vector<double> solve(vector<vector<double>> & D) {
    int m = D.size(), n = D.front().size();

    vector<double> P(n, 0), Q(n, 0), x(n, 0);

    int cmul = 0, cdel = 2;
    P.front() = -D[0][1]/D[0].front();
    Q.front() = D[0].back()/D[0].front();

    for (int i = 1; i < m; ++i) {
        double r = (D[i][i]+D[i][i-1]*P[i-1]);
        ++cmul;
        if (i != m-1) {
            P[i] = -D[i][i+1]/r;
            ++cdel;
        }
        Q[i] = (D[i].back()-D[i][i-1]*Q[i-1])/r;
        ++cmul; ++cdel;
    }

    x.back() = Q.back();
    for (int i = n - 2; i >= 0; --i) {
        x[i] = P[i]*x[i+1]+Q[i];
        ++cmul;
    }

    fout << "cmul = " << cmul << "\tcdel = " << cdel << "\tsum = " << cmul+cdel << endl;

    x.pop_back();
    return x;
}

int main() {
    int n;
    fin >> n;
    vector<vector<double>> s(n, vector<double>(n+1, 0));

    for (auto & row : s)
        for (auto & it : row)
            fin >> it;


    auto a = solve(s);

    fout << "Answer:" << endl;
    for (auto & it : a)
        fout << setw(12) << (fabs(it) > eps ? it : 0);
    fout << endl;

    return 0;
}