Метод простой итерации - путаница в векторе невязки. C++

[Drag0n73]

Вектор невязки в программе увеличивается, хотя по логике должен с каждой итерацией стремиться к нулю. Помогите найти ошибку в коде. Выход из цикла while осуществляется после того, как норма вектора невязки станет меньше заданной константы s

#include <iostream>
#include <iomanip>
#include <cmath>
using namespace std;

int main()
{
	int matrix[4][4]={2,3,-1,-1,5,-1,-1,6,1,7,3,4,-1,-2,1,4};
	int right[4][1]={3,9,15,2};
	int right1[4][1];
	int matrix_t[4][4];
	double matrix_x[4][4];
	double matrix_x1[4];
	double matrix_x2[4]={0,0,0,0};
	double matrix1[4][4]={1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1};
	double matrix_b[4][4];
	double matrix_1[4][4];
	double matrix_g[4];
	double r[4];
	int p=0;
	int it=0;
	double s=0.000000001;
	double m=0;
	double r1=1;
	//транспонируем
	for (int j=0; j<4; j++) {
		for (int i=0;i<4;i++) {
			matrix_t[i][j]=matrix[j][i];
		}
	}
	//умножение матриц
	for (int j=0; j<4; j++) {
		for (int i=0;i<4;i++) {
			p=0;
			for (int k=0; k<4; k++) {
				p+=matrix_t[j][k]*matrix[k][i]; }
			matrix_x[j][i]=p; 
		}
	}
	//умножение матриц
	for (int j=0; j<4; j++) {
		for (int i=0;i<1;i++) {
			p=0;
			for (int k=0; k<4; k++) {
				p+=matrix_t[j][k]*right[k][i]; }
			right1[j][i]=p; 
		}
	}
	//возведение в квадрат
	for (int j=0; j<4; j++) {
		for (int i=0;i<4;i++) {
			m+=matrix_x[j][i]*matrix_x[j][i];
		}
	}
	//2/m
	m=sqrt(m);
	m=2/m;
	//2/m*A
	for (int j=0; j<4; j++) {
		for (int i=0;i<4;i++) {
			matrix_1[j][i]=m*matrix_x[j][i];
		}
	}
	//B
	for (int j=0; j<4; j++) {
		for (int i=0;i<4;i++) {
			matrix_b[j][i]=matrix1[j][i]-matrix_1[j][i];
		}
	}
	//G
	for (int j=0; j<4; j++) {
			matrix_g[j]=m*right1[j][0];
	}
	//начальное приближение
	for (int j=0; j<4; j++) {
			matrix_x1[j]=matrix_g[j];
	}


	while (s<=r1) {
		for (int j=0; j<4; j++) {
				r[j]=right[j][0];
		}
		if (it>0) {
			for (int j=0; j<4; j++) {
				matrix_x1[j]=matrix_x2[j];
			}
		}
		it++;
		for (int j=0; j<4; j++) {
			for (int i=0;i<4;i++) {
				matrix_x2[j]+=matrix_b[i][j]*matrix_x1[i];
			}
			matrix_x2[j]+=matrix_g[j];
		}
		for (int j=0; j<4; j++) {
			for (int i=0;i<4;i++) {
				r[j]=r[j]-(matrix[i][j]*matrix_x2[i]); //вектор невязки
			}
		}
		r1=0;
		for (int i=0;i<4;i++) {
			r1+=r[i]*r[i]; //норма вектора невязки
		}
		r1=sqrt(r1);
	}
	cout<<it<<endl;
	for (int i=0; i<4; i++) {
		cout<<matrix_x2[i]<<" ";
	}
}