Решение СЛАУ методом Гаусса

[doter]

Доброго времени суток . Пишу программу , которая должна решать системы линейных уравнений методом Гаусса(методом ведущего элемента в столбце) . Код написан , но почему-то возникают проблемы "на ровном месте" , например , неправильно находит ведущий элемент , хотя алгоритм вроде бы правильный . Буду благодарен за любую помощь .

#include<iostream.h>
#include<math.h>

void swap_rows(double **, int, int, int);
void swap(double *, int, int);
void div_row(double **, int, int, double);
void Gauss_method(double **, double *, double *, int); 

int main(){
	int i, j, n;
	cout<<"System : Ax=b (A - n x n, b - n x 1)\n";
	cout<<"Enter n : "; 
	cin>>n;
	double ** A = new double * [n];
	for(i = 0; i < n; ++i) {
		A[i] = new double [n];
	}
	double * b = new double[n];
	cout<<"Enter the matrix A : \n";
	for(i = 0; i < n; ++i) {
		for(j = 0; j < n; ++j) {
			cin>>A[i][j];
		}
	}
	cout<<"Enter the vector b : \n";
	for(i = 0; i < n; ++i) {
		cin>>b[i];
	}
	// Решение системы
	double * x = new double [n];
	Gauss_method(A, b, x, n);
	cout<<"The solution : \n";
	for(i = 0; i < n; ++i) {
		cout<<"x_"<<i+1<<" = "<<x[i]<<"\n";
	}
	cout<<endl;
	// Очистка памяти	
	for(i = 0; i < n; ++i) {
		delete  A[i];
	}
    delete [] A;
	delete [] b;
	delete [] x;

	return 0;
}

void swap_rows(double ** A, int n, int i, int j) {
	double tmp = 0;
	for(int k = 0; k < n; ++k) {
		tmp = A[i][k];
		A[i][k] = A[j][k];
		A[j][k] = tmp;
	}
}

void swap(double * b, int i, int j) {
	double tmp = b[i];
	b[i] = b[j];
	b[j] = tmp;
}

void div_row(double ** A, int n, int i, double d) {
	for(int k = 0; k < n; ++k){
		A[i][k] /= d;
	}
}

void Gauss_method(double ** A, double * b, double * x, int n) {
	int i, j, k;
	int row;
	double d = 0, tmp = 0;
	for(k = 0; k < n; ++k) {
		for(i = k; i < n; ++i) {
			if(fabs(A[i][k]) > d){
				d = fabs(A[i][k]);
				row = i;
			}
		}
		if(row != k) {
			swap(b, row, k);
			swap_rows(A, n, row, k);
		}
		b[k] /= A[k][k];
		div_row(A, n, k, A[k][k]);
		for(i = k + 1; i < n; ++i){
			tmp = A[i][k];
			b[i] -= tmp * b[k];
			for(j = k; j < n; ++j) {
				A[i][j] -= tmp * A[k][j];
			}
		}
	}
	for(i = n - 1; i >= 0; --i) {
		x[i] = b[i];
		for(j = i+1; j < n; ++j){
			x[i] -= x[j] * A[i][j];
		}
	}
}