Cod sursa(job #2487634)

Utilizator Kln1000Ciobanu Bogdan Kln1000 Data 5 noiembrie 2019 01:46:10
Problema Algoritmul lui Gauss Scor 10
Compilator cpp-64 Status done
Runda Arhiva educationala Marime 2.21 kb
// Copyright 2019 Ciobanu Bogdan-Calin
// [email protected]

#include <iostream>
#include <fstream>
#include <iomanip>

std::ifstream f("gauss.in");
std::ofstream t("gauss.out");

#define NMAX 310
#define EPS 1e-8

template <typename T>
using Matrix = T[NMAX][NMAX];

template <typename T>
using Vector = T[NMAX];

int main() {
    // n is the number of equations, m the number of unknowns
    Matrix<double> A = {};
    Vector<double> b = {};
    Vector<double> x = {};
    int n, m;

    f >> n >> m;

    auto subtract_line_with_scalar = [&A, &b, m](int dst, int src, double scalar, int start = 0) mutable {
        for (int i = start; i < m; ++i)
            A[dst][i] -= A[src][i] * scalar;

        b[dst] -= b[src] * scalar;
    };

    auto divide_line_with_scalar = [&A, &b, m](int dst, double scalar, int start = 0) mutable {
        for (int i = start; i < m; ++i)
            A[dst][i] /= scalar;

        b[dst] /= scalar;
    };

    auto not_linear_independent = [&A, m](int dst) {
        bool flag = true;
        for (int i = 0; i < m; ++i) {
            std::cout << A[dst][i] << " ";
            flag &= (A[dst][i] * A[dst][i] < EPS ? true : false);
        }
        std::cout << std::endl << flag << std::endl;
        return flag;
    };

    // reading the matrix and results
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < m; ++j)
            f >> A[i][j];
        f >> b[i];
    }

    for (int i = 0; i < m; ++i) {  // for each column
        divide_line_with_scalar(i, A[i][i], i); // no pivoting, but scaling
        for (int j = i + 1; j < n; ++j) { // rest to 0
            if (not_linear_independent(j)) {
                t << "Imposibil\n";
                return 0;
            }
            subtract_line_with_scalar(j, i, A[j][i], i);
        }

    }

    double sigma;
    // now perform back substitution
    for (int i = m - 1; i >= 0; --i) {
        sigma = .0;
        for (int j = m - 1; j > i; --j) {
            sigma += x[j] * A[i][j];
        }
        x[i] = (b[i] - sigma) / A[i][i];
    }

    t << std::setprecision(10) << std::fixed;
    for (int i = 0; i < m; ++i)
        t << x[i] << " ";

    return 0;
}