【数值分析练习】三阶矩阵jacobi迭代法

其实是纯c语言，多封装几个函数即可。 

//n*n矩阵和n*1矩阵的jacobi迭代法
#include
#include
#define ROW 3
#define COL 3
void input(double A[][COL + 1], double B[][2], double G[][COL + 1], double F[][2])
{
	//读取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]);
	//逐行对G，F赋值
	double a_ii;
	for (int i = 1; i <= ROW; i++)
	{
		//G
		a_ii = A[i][i];
		if (a_ii == 0) //错误情况
		{
			puts("a_ii==0");
			return;
		}
		for (int j = 1; j <= COL; j++)
		{
			if (j == i) //对角线0
				G[i][j] = 0;
			else //其他为正常
				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("%ft", M[i][j]);
		putchar('n');
	}
}
void printVector(double V[][2])
{
	for (int i = 1; i <= ROW; i++)
		printf("%fn", 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];//采取从1开始记法
	double GX[ROW + 1][2], AX[ROW + 1][2];
	input(A, B, G, F);

	//迭代过程
	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);//作差分析误差
		printVector(X);
		putchar('n');
	} while (getVectorNorm(R) > 0.00001);//范数精度

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

	return 0;
}