[数值分析]线性方程组求解：列主元高斯消元法

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

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

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

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

	//计算线性方程组的解
	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);
	MatrixPrint(Mat_Solve, n, 1);

	system("Pause");
	return 0;
}

double* MatrixSolve(double* arr_co_in, double* arr_y_in, const int n)
{
	int arr_co_len = n * n;
	int i;
	int col_scan, row_scan, resi_row_scan;
	double pivot, y_row_scan;//主元

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

	//复制一份拷贝以防止改变原有数组
	double *arr_co = (double *)malloc(arr_co_len * sizeof(double));
	if (arr_co == NULL) return NULL;
	for (i = 0; i < arr_co_len; i++) arr_co[i] = arr_co_in[i];

	double *arr_y = (double *)malloc(n * sizeof(double));
	if (arr_y == NULL) return NULL;
	for (i = 0; i < n; i++) arr_y[i] = arr_y_in[i];

	//开始计算
	for (row_scan = 0; row_scan < n; row_scan++)
	{
		//列主元换位
		MatrixPivotExchange(arr_co, arr_y, n, row_scan);

		//原矩阵与增广矩阵行元素均需要除以列主元
		pivot = arr_co[row_scan * (n + 1)];

		for (col_scan = row_scan; col_scan < n; col_scan++) // 下三角
		{
			arr_co[row_scan * n + col_scan] = arr_co[row_scan * n + col_scan] / pivot;
		}
		arr_y[row_scan] = arr_y[row_scan] / pivot;

		//通过主元行消除其余行首位元素
		for (resi_row_scan = row_scan + 1; resi_row_scan < n; resi_row_scan++)
		{
			pivot = arr_co[resi_row_scan * n + row_scan];

			for (col_scan = row_scan; col_scan < n; col_scan++) //下三角
			{
				arr_co[resi_row_scan * n + col_scan] = arr_co[resi_row_scan * n + col_scan] - pivot * arr_co[row_scan * n + col_scan];
			}

			arr_y[resi_row_scan] = arr_y[resi_row_scan] - pivot * arr_y[row_scan];
		}
	}

	//逐步通过三角矩阵求解
	for (row_scan = n - 1; row_scan >= 0; row_scan--)
	{
		y_row_scan = arr_y[row_scan];

		for (col_scan = row_scan + 1; col_scan < n; col_scan++)
		{
			y_row_scan = y_row_scan - arr_co[row_scan*n + col_scan] * Mat_Solve[col_scan];
		}
		Mat_Solve[row_scan] = y_row_scan / arr_co[row_scan * (n + 1)];
	}

	free(arr_co);
	free(arr_y);
	return Mat_Solve;
}

//列主元换位
void MatrixPivotExchange(double* arr, double* arr_y, const int n, const int row) //row为现在处理的行
{
	int    row_scan, col_scan;						//循环标记
	int	   col				= row;					//列主元的扫描列 就是输入的行
	double max_value		= arr[row * (n + 1)];	//输入时行的列主元值
	int	   max_value_row	= row;					//记录最大列主元的所在行
	double cmp_calue		= 0;
	double temp				= 0;

	int org_index, max_index;

	//找到列主元的最大值所在行
	for (row_scan = row + 1; row_scan < n; row_scan++)
	{
		cmp_calue = arr[row_scan * n + col];
		if (max_value < cmp_calue)
		{
			max_value	  = cmp_calue;
			max_value_row = row_scan;
		}
	}

	//原矩阵与逆阵均需要换行
	if (row != max_value_row)
	{
		for (col_scan = 0; col_scan < n; col_scan++)
		{
			org_index = row * n + col_scan;
			max_index = max_value_row * n + col_scan;

			temp = arr[org_index];
			arr[org_index] = arr[max_index];
			arr[max_index] = temp;

		}

		//兼容
		if (arr_y != NULL)
		{
			temp = arr_y[row];
			arr_y[row] = arr_y[max_value_row];
			arr_y[max_value_row] = temp;
		}
	}
	
	return;
}


//矩阵打印方法
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;
}