列主元消去.cpp

求线性方程组的解和系数矩阵的逆

#include <stdio.h>
#include <math.h>
#include <process.h>

void Gaussian_elimination(float * x, float * coefficient, int width, int height);
int max(float * coefficient, int width, int height, int i);
void exchange(float * coefficient, int width, int height, int order, int i);
void inverse_matrix(float * coefficient, float * inverse_coefficient, int dimension);

void main()
{
	int width, height, dimension;
	int i, j;
    float * coefficient, * inverse_coefficient, * x;

	printf("/******************************************************************************/\n");
	printf("                                   程序说明\n");
	printf("\n 该程序用于解齐次线性方程组的解，输入为增广矩阵，输出为解向量和系数矩阵的逆矩阵\n");
	printf("\n/******************************************************************************/\n");
	printf("\n输入方程的维数:\n");
	scanf("%d", &dimension);
    
	height = dimension;
	width = dimension + 1;

	coefficient = new float[width * height];
	inverse_coefficient = new float [dimension * dimension];
    x = new float[height];

	for (i = 0; i < height; i++)
	{
		printf("输入矩阵第%d行:\n", i + 1);
		for (j = 0; j < width; j++)
	        scanf("%f", &coefficient[i * width + j]);
	}
    //用高斯消去法求方程的解
	Gaussian_elimination(x, coefficient, width, height);
    //用高斯-约旦消去法求矩阵的逆
    inverse_matrix(coefficient, inverse_coefficient, dimension);

	printf("方程的解为:\n");
	for (i = 0; i < dimension; i++)
	    printf("x%d = %f\n", i + 1, x[i]);

	printf("系数矩阵的逆矩阵为:\n");
	for (i = 0; i < dimension; i++)
	{
		for (j = 0; j < dimension; j++)
			printf("%f ", inverse_coefficient[i * dimension + j]);
		printf("\n");
	}
	system("PAUSE");

	delete[] coefficient;
	delete[] inverse_coefficient;
}
//***************************************************************************************
//函数功能: 用高斯列主元消去法求方程的解
//函数输入: 增广矩阵的系数 * coefficient, 增广矩阵的高度 height, 以及增广矩阵的宽度 width
//函数输出: 解向量 x
//ps:要求增广矩阵的高度 height就是其所代表方程的维数,而矩阵的宽度就是width = height + 1
//***************************************************************************************
void Gaussian_elimination(float * x, float * coefficient, int width, int height)
{
	int i, j, k;
    int order;
	float aii, rate, temp;
    float *temp_matrix = new float [width * height];

	for (i = 0; i < height; i++)
		for (j = 0; j < width; j++)
			temp_matrix[i * width + j] = coefficient[i * width + j];

    for (i = 0; i < height - 1; i++)
	{
		order = max(temp_matrix, width, height, i);
		exchange(temp_matrix, width, height, order, i);
        aii = temp_matrix[i * width + i];

		for (j = i + 1; j < height; j++)
		{
			rate = - temp_matrix[j * width + i] / aii;

			for (k = i; k < width; k++)
			{
				temp_matrix[j * width + k] = temp_matrix[i * width + k] * rate 
					                       + temp_matrix[j * width + k];
			}
		}
	}
    
	//回代,先算最后一个 x
	x[height - 1] = temp_matrix[height * width - 1] / temp_matrix[height * width - 2];

	//iteration
	for (i = height - 2; i >= 0; i--)
	{
		// ∑akl * xl
		temp = 0;
		for (j = i + 1; j < width - 1; j++)
            temp += temp_matrix[i * width + j] * x[j];

		x[i] = (temp_matrix[(i + 1) * width - 1] - temp) / temp_matrix[i * width + i];
	}

	delete [] temp_matrix;
}
//**************************************************************************************
//函数功能: 求第i列中绝对值最大的数所在的行数
//函数输入: 增广矩阵的系数 coefficient, 矩阵的宽度 width, 高度 height, 以及所要求的所在列i
//函数输出: 最大值所在的行数
//ps:要求增广矩阵的高度 height就是其所代表方程的维数,而矩阵的宽度就是width = height + 1
//**************************************************************************************
int max(float * coefficient, int width, int height, int i)
{
	float temp;
	int j, ret;

	temp = coefficient[i * width + i];
    ret = i;

	for (j = i; j < height; j++)
	{
		if (fabs(coefficient[j * width + i]) > fabs(temp))
		{
			temp = coefficient[j * width + i];
			ret = j;
		}
	}

	return ret;
}
//****************************************************************************************
//函数功能: 将矩阵的两行交换
//函数输入: 增广矩阵的系数 coefficient, 矩阵的宽度 width, 高度 height, 以及所要交换的行 i和 order
//函数输出: 交换后的系数矩阵
//****************************************************************************************
void exchange(float * coefficient, int width, int height, int order, int i)
{
	float temp;
	int j;

	for (j = i; j < width; j++)
	{
		temp = coefficient[i * width + j];
        coefficient[i * width + j] = coefficient[order * width + j];
        coefficient[order * width + j] = temp;
	}
}
//******************************************************************
//函数功能: 用列主元高斯-约旦消去法求矩阵的逆
//函数输入: 增广矩阵的系数 coefficient, 矩阵的维数 dimension
//函数输出: 系数矩阵的逆矩阵 inverse_coefficient
//ps: 对于输入是一般方阵的情况,下面的注释说明了如何更改程序
//******************************************************************
void inverse_matrix(float * coefficient, float * inverse_coefficient, int dimension)
{
	int i, j, k;
    int order;
	float aii, rate;
    float * temp_matrix = new float [dimension * dimension * 2];
	
	//initial
	for (i = 0; i < dimension; i++)
		for (j = 0; j < dimension * 2; j++)
			if(j < dimension)//这是针对的增广矩阵的系数矩阵求逆的,所以后面的是coefficient[i * (dimension + 1) + j]
				//若输入是一般的方阵,则应变成coefficient[i * dimension + j]
				temp_matrix[i * dimension * 2 + j] = coefficient[i * (dimension + 1) + j];
			else
			{
				if ((j - dimension) == i)
				    temp_matrix[i * dimension * 2 + j] = 1;
				else
                    temp_matrix[i * dimension * 2 + j] = 0;
			}

    for (i = 0; i < dimension; i++)
	{
		order = max(temp_matrix, dimension * 2, dimension, i);
		exchange(temp_matrix, dimension * 2, dimension, order, i);
        aii = temp_matrix[i * dimension * 2 + i];

		for (j = 0; j < dimension; j++)
		{
			if (j != i)
			{
				rate = - temp_matrix[j * dimension * 2 + i] / aii;

				for (k = i; k < dimension * 2; k++)
				{
					temp_matrix[j * dimension * 2 + k] = temp_matrix[i * dimension * 2 + k] * rate 
											   + temp_matrix[j * dimension * 2 + k];
				}
			}
		}
	}

    for (i = 0; i < dimension; i++)
		for (j = 0; j < dimension; j++)
			inverse_coefficient[i * dimension + j] = temp_matrix[i * dimension * 2 + j + dimension] 
			                                       / temp_matrix[i * dimension * 2 + i];

	delete []temp_matrix;
}