[数值分析]线性方程组求解：Jacobi迭代法

// ConsoleApplication1.cpp : 定义控制台应用程序的入口点。
//

#include "stdafx.h"
#include <stdlib.h> 
#include "stdio.h"
#include "math.h"

void MatrixPrint(double* arr, const int row, const int col);
double* MatrixSolve(double* arr_co_in, double* arr_y_in, const int n);

int _tmain(int argc, _TCHAR* argv[])
{   
	int i,p;
	//统一使用一维数组实现
	//Matlab 解为 [6 -8 11]
	const int n = 3;
	double Mat_A[n * n] = { 5, 2, -2, 1, 3, 1, 2, 2, 5 }; //矩阵
	double Mat_Y[n]		= { 3, 6, 3 };					 //增广矩阵

	//计算线性方程组的解
	double *Mat_Solve = MatrixSolve(Mat_A, Mat_Y, n);

	//打印增广矩阵 和 方程组的解
	double *Mat_print = (double *)malloc((n*(n+1)) * sizeof(double));
	for (i = 0, p = 0; i < n * n; i++)
	{	
		Mat_print[p++] = Mat_A[i];
		if((i+1) % 3 == 0)
		{
			Mat_print[p] = Mat_Y[(int)(((i+1) / 3) - 1)];
			p++;
		}
	}

	MatrixPrint(Mat_print, n, n + 1);
	if (Mat_Solve != NULL)
		MatrixPrint(Mat_Solve, n, 1);
	else
		printf("矩阵不收敛");

	system("Pause");
	return 0;
}

double* MatrixSolve(double* arr_co_in, double* arr_y_in, const int n)
{
	const double iter_presion = 0.00001; //迭代精度
	const int    iter_t_limit = 100;     //迭代最长次数
	int iter_t = 0;						 //迭代次数
	int arr_len = n * n;
	int i, j, count;
	double sum_temp, max_temp;

	//通过无穷范数（谱半径）计算是否会收敛
	//申请缓存

	for (i = 0, max_temp = 0; i < n; i++)
	{
		sum_temp = 0;
		for (j = 0; j < n; j++)
		{	
			if (i != j)
			{
				sum_temp = sum_temp + fabs(arr_co_in[i * n + j] / arr_co_in[i * (n + 1)]);
			}
		}
		if (sum_temp > max_temp)
			max_temp = sum_temp;
	}
	if (max_temp >= 1) return NULL;

	//申请解的缓存
	double *Mat_Solve = (double *)malloc(n * sizeof(double));
	if (Mat_Solve == NULL) return NULL;

	//前解用于检验收敛
	double *Mat_Solve_Pre = (double *)malloc(n * sizeof(double));
	if (Mat_Solve_Pre == NULL) return NULL;
	
	//初值均赋值为0
	for (i = 0; i < n; i++)
	{
		Mat_Solve[i] = 0;
		Mat_Solve_Pre[i] = 0;
	}

	while (1)
	{
		//迭代限制
		iter_t++;
		if (iter_t > iter_t_limit)
		{
			return NULL;
		}

		//进行一次迭代
		for (i = 0; i < n; i++)
		{
			sum_temp = arr_y_in[i];
			//每条线性方程组移向右侧并除以主元系数
			for (j = 0; j < n; j++)
			{
				if (i != j)
				{
					//系数全部使用旧值
					sum_temp = sum_temp - Mat_Solve_Pre[j] * arr_co_in[i * n + j];
				}
			}
			sum_temp = sum_temp / arr_co_in[i * (n + 1)];
			Mat_Solve[i] = sum_temp;
		}

		//检验迭代收敛
		for (i = 0, count = 0; i < n; i++)
		{
			if (fabs(Mat_Solve[i] - Mat_Solve_Pre[i]) <= iter_presion)
			{
				count++;
			}
		}

		if (count == n) 
			break;
		else
		{
			for (i = 0, count = 0; i < n; i++)
			{
				Mat_Solve_Pre[i] = Mat_Solve[i];
			}
		}
	}

	free(Mat_Solve_Pre);
	return Mat_Solve;
}



//矩阵打印方法
void MatrixPrint(double* arr, const int row, const int col)
{
	int len = row * col;
	int i,col_count;

	for (i = 0, col_count = 0; i < len; i++)
	{
		printf("%f\t", arr[i]);

		//单换行
		if (++col_count >= col)
		{
			printf("\n");
			col_count = 0;
		}
	}

	//跳空换行
	printf("\n");

	return;
}