// метод Гаусса-Жордана с выбором ведущего элемента по столбцу
// a - матрица, b - правый столбец, потом решение
void gauss(double *a, double *b, int n) {
    int i, j, k, mxi;
    double mx,w;
    for (i = 0; i < n; i++) 
    {
        mx = fabs(a[i*n+i]);
        mxi = i;
        for (k = i + 1; k < n; k++)
          if ( fabs(a[k*n+i]) > mx ) { mx = fabs(a[k*n+i]); mxi = k; }
// меняем местами строки:
        for (j = i; j < n; j++) {w=a[i*n+j]; a[i*n+j]=a[mxi*n+j]; a[mxi*n+j]=w;}
        w=b[i]; b[i]=b[mxi]; b[mxi]=w; 
// нормализация строки:        
        for (j = i+1; j<n; j++) a[i*n+j]/=a[i*n+i];
        b[i]/=a[i*n+i]; 
        a[i*n+i]=1.0;
        
        for (k = 0; k < n; k++) {
			if(k != i){
				for (j = i + 1; j < n; j++) a[k*n+j] -= a[k*n+i]*a[i*n+j];
                b[k] -= a[k*n+i]*b[i];
				a[k*n+i] = 0.0;
			}
        }
    }
}

// a - матрица, x - решение, b - исходная правая часть, db=b-Ax - невязка
void delta_b(double *a, double *x, double *b, double *db, int n) 
{
	int i,j;	
	for(i=0;i<n;i++){
		db[i]=b[i];
		for(j=0;j<n;j++) db[i]-=a[i*n+j]*x[j];
	}
}