Cod sursa(job #1283170)

Utilizator ucc_5Usurelu Catalin ucc_5 Data 5 decembrie 2014 10:55:13
Problema Algoritmul lui Gauss Scor 0
Compilator cpp Status done
Runda Arhiva educationala Marime 2.94 kb
#include <fstream>
#include <iostream>
#include <iomanip> 
#include <algorithm>
#include <cmath>

#define EPS 0.0000001

using namespace std;

ofstream g("gauss.out");

void debug(double**& A, const int& N, const int& M)
{
	for(int i = 1; i <= N; i++, cout << endl)
        for(int j = 1; j <= M + 1; j++)
        	cout << A[i][j] << " ";

    cout << endl;
}

void read_input(double**& A, int& N, int& M)
{
	ifstream f("gauss.in");
	f >> N >> M;

	A = new double*[N + 1];
	for(int i = 1; i <= N; i++)
	{
		A[i] = new double[M + 2];
	}

	for(int i = 1; i <= N; i++)
	{
		for(int j = 1; j <= M + 1; j++)
		{
			f >> A[i][j];
		}
	}

	f.close();
}

int argmax(double**& A, const int& N, const int& col)
{
	int row_max = col;
	for(int j = col + 1; j <= N; j++)
	{
		if(abs(A[j][col]) > abs(A[row_max][col]))
		{
			row_max = col;
		}
	}

	return row_max;
}

void swap_rows(double**& A, const int& N, const int& M, const int& k, const int& pivot_index)
{
	for(int i = 1; i <= M + 1; i++)
	{
		swap(A[k][i], A[pivot_index][i]);
	}
}

// Compute in-place row echelon form
bool compute_row_echelon_form(double**& A, const int& N, const int& M)
{
	//debug(A, N, M);

	for(int k = 1; k < min(N, M); k++)
	{
		int pivot_index = argmax(A, N, k);

		//cout << "Pivot index " << pivot_index << endl;

		if(A[pivot_index][k] > -EPS && A[pivot_index][k] < EPS)
		{
			return false;
		}

		if(k != pivot_index)
		{
			swap_rows(A, N, M, k, pivot_index);
		}

		for(int i = k + 1; i <= N; i++)
		{
			for(int j = k + 1; j <= M + 1; j++)
			{
				A[i][j] = A[i][j] - A[k][j] * (A[i][k] / A[k][k]);
			}
			A[i][k] = 0;
		}

		//debug(A, N, M);
	}

	return true;
}

bool compute_unknowns(double*& X, double**& A, const int& N, const int& M)
{
	//cout << "compute_unknowns " << endl;

	 //Calculul necunoscutelor.
    for(int i = N; i >= 1; i--)
    {
        for(int j = 1; j <= M + 1; j++)
        {
            if(A[i][j] < -EPS || A[i][j] > EPS)
            {
                // Toti coeficientii sunt 0 iar rezultatul este diferit de 0.
                if(j == M + 1) // adica ultima coloana, b-ul
                {
                    return false;
                }
 
                // Calculam pe necunoscuta j = rezultatul ecuatiei
                // i - necunoscutele j+1...M inmultite cu coeficientii lor din aceasta ecuatie.
                X[j] = A[i][M + 1];
                for(int k = j + 1; k <= M; k++)
                {
                    X[j] -= X[k] * A[i][k];
                }

                X[j] /= A[i][i];
 
                break;
            }
        }
    }

    return true;
}


int main()
{
	int N, M;
	double** A;
	double* X;

	read_input(A, N, M);

	X = new double[M + 2];

	bool ret = compute_row_echelon_form(A, N, M);

	if(!ret) 
	{
		cout << "Imposibil\n";
	}
	else
	{
		ret = compute_unknowns(X, A, N, M);

		if(!ret) 
		{
			cout << "Imposibil\n";
		}
		else
		{
			for(int i = 1; i <= M; i++)
			{
				cout << std::fixed << std::setprecision(10) << X[i] << " ";
			}
		}
	}

	return 0;
}