[numerical analysis exercise] Jacobi iteration method of third-order matrix

Actually, it's pure c Language , Just encapsulate a few more functions . 

//n*n Matrix and n*1 Matrix jacobi Iterative method 
#include<cstdio>
#include<cmath>
#define ROW 3
#define COL 3
void input(double A[][COL + 1], double B[][2], double G[][COL + 1], double F[][2])
{
	// Read A,B
	for (int i = 1; i <= ROW; i++)
		for (int j = 1; j <= COL; j++)
			scanf("%lf", &A[i][j]);
	for (int i = 1; i <= ROW; i++)
		scanf("%lf", &B[i][1]);
	// Line by line G,F assignment 
	double a_ii;
	for (int i = 1; i <= ROW; i++)
	{
		//G
		a_ii = A[i][i];
		if (a_ii == 0) // Wrong situation 
		{
			puts("a_ii==0");
			return;
		}
		for (int j = 1; j <= COL; j++)
		{
			if (j == i) // Diagonals 0
				G[i][j] = 0;
			else // Others are normal 
				G[i][j] = -A[i][j] / a_ii;
		}
		//F
		F[i][1] = B[i][1] / a_ii;
	}

}
void printMatrix(double M[][COL + 1])
{
	for (int i = 1; i <= ROW; i++)
	{
		for (int j = 1; j <= COL; j++)
			printf("%f\t", M[i][j]);
		putchar('\n');
	}
}
void printVector(double V[][2])
{
	for (int i = 1; i <= ROW; i++)
		printf("%f\n", V[i][1]);
}
double getVectorNorm(double V[][2])
{
	double max = 0, temp;
	for (int i = 1; i <= ROW; i++)
	{
		temp = fabs(V[i][1]);
		if (temp > max)
			max = temp;
	}
	return max;
		
}
void vectorMinus(double R[][2], double Va[][2], double Vb[][2])
{
	for (int i = 1; i <= ROW; i++)
		R[i][1] = Va[i][1] - Vb[i][1];
}
void matrixMult(double R[][2], double G[][4], double X[][2])
{
	for (int i = 1; i <= ROW; i++)
	{
		for (int j = 1; j <= 1; j++)
		{
			double count = 0;
			for (int k = 1; k <= COL; k++)
				count += G[i][k] * X[k][j];
			R[i][j] = count;
		}
		
	}
}
void vectorAdd(double R[][2], double GX[][2], double F[][2])
{
	for (int i = 1; i <= ROW; i++)
		R[i][1] = GX[i][1] + F[i][1];
}
void vectorCopy(double A[][2], double C[][2])
{
	for (int i = 1; i <= ROW; i++)
		C[i][1] = A[i][1];
}
int main(void)
{
	freopen("input.txt", "r", stdin);
	double A[ROW + 1][COL + 1], B[ROW + 1][2];
	double G[ROW + 1][COL + 1], F[ROW + 1][2];
	double pre_X[ROW + 1][2], X[ROW + 1][2],R[ROW + 1][2];// Take from 1 Start remembering 
	double GX[ROW + 1][2], AX[ROW + 1][2];
	input(A, B, G, F);

	// iterative process 
	vectorCopy(F, X);
	vectorCopy(F, pre_X);
	do {
		vectorCopy(X, pre_X);
		matrixMult(GX, G, X);
		vectorAdd(X, GX, F);
		vectorMinus(R, X, pre_X);// Error of error analysis 
		printVector(X);
		putchar('\n');
	} while (getVectorNorm(R) > 0.00001);// Norm accuracy 

	puts("last result X:");
	printVector(X);
	puts("matching  result AX:");
	matrixMult(AX, A, X);
	printVector(AX);

	return 0;
}