#include stdafx.h

#include "stdafx.h"

#include <iostream>

#include <iomanip>

#include <cmath>

#include <iostream>

#include <iomanip>

#include <cmath>

using namespace std;

void enterMatrix (double**, int n);

void enterVector (double*, int n);

void printMatrix (double**, int n);

void printVector (double*, int n);

void solveGaussMethod (double** A, double* B, double* x, int n, bool& isSolved);

int indexMax(double** A, int k, int n);

int main()

{

cout << "Enter please system dimension: ";

int n = 0;

cin >> n;

double** a = new double* [n];

for (int i = 0; i < n; i++)

a[i] = new double [n];

double* b = new double[n];

double* x = new double[n];

enterMatrix(a, n);

enterVector(b, n);

printMatrix(a, n);

printVector(b, n);

int k=0;

bool isSolved = true;

solveGaussMethod(a, b, x, n, isSolved);

for (int i = 0; i < n; i++)

delete [] a[i];

delete [] a;

delete [] b;

delete [] x;

return 0;

}

void enterMatrix (double** a, int n)

{

cout << "Enter please matrix" << endl;

for (int i = 0; i < n; i++)

for (int j = 0; j < n; j++)

{

cout << "A[" << i << "][" << j << "] = ";

cin >> a[i][j];

}

}

void enterVector (double* b, int n)

{

cout << "Enter please vector" << endl;

for (int i = 0; i < n; i++)

{

cout << "B[" << i << "] = ";

cin >> b[i];

}

}

void printMatrix (double** a, int n)

{

cout << "\n ******* Matrix *****" << endl;

for (int i = 0; i < n; i++)

{

for (int j = 0; j < n; j++)

cout << setw(5) << a[i][j];

cout << endl;

}

}

void printVector (double*b, int n)

{

cout << "***** Vector *******" << endl;

for (int i = 0; i < n; i++)

cout << setw(12) << b[i];

cout<<endl;

}

void solveGaussMethod (double** A, double* B, double* x, int n, bool& isSolved)

{

double** a = new double* [n];

for (int i = 0; i < n; i++)

a[i] = new double [n];

for (int i = 0; i < n; i++)

for (int j = 0; j < n; j++)

a[i][j] = A[i][j];

double* b = new double[n];

for (int j = 0; j < n; j++)

b[j] = B[j];

// Прямой ход

isSolved = true;

for(int k = 0; k < n; k++)

{

//cout << "indMax = " << indexMax(A, k, n) << endl;

//printMatrix(A, n);

//printVector(B, n);

int indMax = indexMax(A, k, n);

double* temp = A[k];

A[k] = A[indMax];

A[indMax] = temp;

double t = B[k];

B[k] = B[indMax];

B[indMax] = t;

//printMatrix(A, n);

//printVector(B, n);

double MAIN = A[k][k];

if (fabs(MAIN) < 1e-15)

{

cout << "Error 1!" << endl;

isSolved = false;

return;

}

for(int j = k; j < n; j++)

A[k][j] = A[k][j]/MAIN;

B[k] = B[k] / MAIN;

for (int i = k+1; i < n; i++)

{

double MAINi = A[i][k];

for (int j = k; j < n; j++)

A[i][j] = A[i][j] - MAINi*A[k][j];

B[i] = B[i] - MAINi*B[k];

}

printMatrix(A, n);

printVector(B, n);

}

// Обратный ход

x[n-1] = B[n-1];

for(int i = n-2; i >= 0; i--)

{

double sum = 0.0;

for(int j = i + 1; j < n; j++)

sum += A[i][j] * x[j];

x[i] = B[i] - sum;

}

printVector(x, n);

// Вектор невязки для Ax-B=F

double summ = 0.0;

for(int j = 0; j < n; j++)

summ += A[0][j] * x[j];

double maxf = fabs(summ - B[0]);

for(int i = 0; i < n; i++)

{

double sum = 0.0;

for(int j = 0; j < n; j++)

sum += A[i][j] * x[j];

double f = sum - B[i];

if (fabs(f) > maxf)

maxf = fabs(f);

cout<< "f = "<< f <<endl;

}

cout<< "norm = "<< maxf <<endl;

}

int indexMax(double** A, int k, int n)

{

double max = fabs(A[k][k]);

int indMax = k;

for (int i = k + 1; i < n; i++)

{

if (max < fabs(A[i][k]))

{

max = fabs(A[i][k]);

indMax = i;

}

}

return indMax;

}

Дата добавления: 2015-11-05