高斯消元法求逆矩阵

  有多组测试数据。每组测试数据先输入一个整数n,表示方阵的阶。然后下面输入n阶方阵。输出其逆矩阵。若无逆矩阵,则输出No inverse matrix。

#include <iostream>
#include <cmath>
#include <algorithm>

using namespace std;

const double eps = 1e-6;

bool is_zero( const double num )
{
	return fabs(num) < eps;
}

void create( double ** & matrix, const int n )
{
	matrix = new double* [n];
	for ( int i = 0; i < n; ++i )
		matrix[i] = new double[n];
}

void input ( double ** matrix, const int n )
{
	for ( int i = 0; i < n; ++i )
	{
		for ( int j  = 0; j < n; ++ j )
			cin >> matrix[i][j];
	}
}

bool inverse ( double ** matrix1, double ** matrix2, const int n )
{
	int i, j;
	for ( i = 0; i < n; ++ i )
	{
		for ( j = 0; j < n; ++ j )
		{
			if ( i == j )
				matrix2[i][j] = 1;
			else
				matrix2[i][j] = 0;
		}
	}
	for ( i = 0; i < n; ++i )
	{
		int rowmaxpos = i;
		for ( j = i + 1; j < n; ++j )
		{
			if ( matrix1[i][j] > matrix1[i][rowmaxpos] )
				rowmaxpos = j;
		}
		for ( j = i; j < n; ++ j )
		{
			swap( matrix1[j][rowmaxpos], matrix1[j][i]);
			swap( matrix2[j][rowmaxpos], matrix2[j][i]);
		}
		if ( !is_zero(matrix1[i][i]) )
		{
			int divisor = matrix1[i][i];
			for ( j = i; j < n; ++ j )
			{
				matrix1[i][j] /= divisor;
				matrix2[i][j] /= divisor;
			}
			for ( j = i + 1; j < n; ++ j )
			{
				int multiple = matrix1[j][i];
				for ( int k = i; k < n; ++ k )
				{
					matrix1[i][j] -= matrix1[i][k] * multiple;
					matrix2[i][j] -= matrix2[i][k] * multiple;
				}
			}
		}
		else
			return false;
	}
	return true;
}

void output( double ** matrix, const int n )
{
	for ( int i = 0; i < n; ++i )
	{
		for ( int j = 0; j < n; ++ j )
			cout << matrix[i][j] << ' ';
		cout<<endl;
	}
}

void destroy( double ** matrix, const int n )
{
	for ( int i = 0; i < n; ++ i )
		delete [] matrix[i];
	delete [] matrix;
}

int main()
{
	int n;
	double ** matrix1;
	double ** matrix2;
	while ( cin >> n )
	{
		create( matrix1, n );
		create( matrix2, n );
		input( matrix1, n);
		if ( inverse(matrix1, matrix2, n) )
			output( matrix2, n );
		else
			cout << "No inverse matrix" << endl;
		destroy( matrix1, n );
		destroy( matrix2, n );
	}
	return 0;
}

点赞