Написал код реализующий метод Гаусса с полным выбором ведущего элемента, но ответ выдает не правильный. Помогите найти ошибку
Код:
#include "stdafx.h"
#include <conio.h>
#include <stdio.h>
#include <iostream>
#include <iomanip>
#include <math.h>
#define N 6
void PrintArray(double A[][N]){ //вывод матрицы
for(int i=0; i<N; i++){
for(int j=0; j<N; j++)
std::cout<< A[i][j] <<" ";
std::cout<<std::endl;
}
}
void PrintVector(double X[]){ // вывод вектора
for(int i=0; i<N; i++){
std::cout << std::fixed << std::setw(22) << std::setprecision(20) << X[i];
std::cout<<std::endl;
}
}
void GilbertM(double A[][N]){ //матрица Гильберта
int i;
int j;
for (i = 0;i<N;i++){
for(j=0;j<N;j++){
A[i][j]=1./(i+1+j);
}
}
}
void VectorB(double B[]){ // заполняем вектор B
int i;
B[0] = 1;
for (i = 1;i<N;i++){
if (i % 2 == 0){
B[i]=-(B[i-1]-1);
}
else{ B[i] = -(B[i - 1] + 1);}
}
}
void SwapColumns(double A[][N],int c1,int c2){ // меняем столбцы матрицы
double temp;
int j;
for(j=0; j<N; j++)
{
temp=A[j][c1];
A[j][c1]=A[j][c2];
A[j][c2]=temp;
}
}
void SwapLines(double A[][N], double B[],int l1,int l2){ // меняем строки матрицы и вектора
double temp;
int j;
for(j=0; j<N; j++)
{
temp=A[l1][j];
A[l1][j]=A[l2][j];
A[l2][j]=temp;
temp=B[l1];
B[l1]=B[l2];
B[l2]=temp;
}
}
void Destroy(double A[][N], double B[], int d){ // зануляем элементы под ведущим элементом
int i,k;
double v;
for(i=d+1; i<N; i++){
v=(-A[i][d]/A[d][d]);
for(k=d; k<N; k++){
A[i][k]=A[i][k]+v*A[d][k];
}
B[i]=B[i]+v*B[d];
}
}
void Solution(double A[][N], double B[], double X[]){ //собственно само решение
int i,j,s,m,k,l,temp;
double max;
int memory[N];
double Y[N];
double v;
for (i=0; i<N; i++)
memory[i]=i;
for (i=0; i<N; i++){
max=0;
k=i;
l=i;
for (s = i;s<N;s++){//ищем максимум в подматрице
for (m = i;m<N;m++){
if (fabs(A[s][m])>max){max=fabs(A[s][m]); k=s; l=m;}
}
}
// if (k!=i)
SwapLines(A,B,i,k);//меням строки
// if (l!=i){
SwapColumns(A,i,l);// меням столбцы и запоминаем какие переменные поменялись
temp = memory[i];
memory[i] = memory[l];
memory[l] = temp;
// }
Destroy(A,B,i); // зануление
}
for (i=N-1; i>-1; i--){ //обратная подстановку
v=0;
for (j=N-1; j>i; j--)
v=A[i][j]*Y[j]+v;
Y[i]=(B[i]-v)/A[i][i];
}
for (i = 0; i < N; i++)//меняем местами переменные как должно быть
X[memory[i]] = Y[i];
}
void F(double A[][N], double X[N], double B[N]){// подставляем A,B,X в уравнение AX-B и смотрим чему равен максимальный элемент
double L[N];
double max;
int i,j;
max=0;
for(i=0; i<N; i++){
L[i]=0;
for(j=0; j<N; j++){
L[i]=L[i]+A[i][j]*X[j];
}
L[i]=L[i]-B[i];
}
for (i=0; i<N; i++){
max=max+fabs(L[i]);
}
std::cout << std::fixed << std::setw(22) << std::setprecision(20) << max<<std::endl<<std::endl;
}
void main(){
double B[N];
double X[N];
double A[N][N];
double max;
GilbertM(A);
VectorB(B);
Solution(A,B,X);
PrintVector(X);
GilbertM(A);
VectorB(B);
std::cout<<std::endl;
F(A,X,B);
getch();
}
