矩阵求逆.cpp

来自「数值计算的例子和参考书籍」· C++ 代码 · 共 47 行

double a[MAX][MAX], I[MAX][MAX];
int n;
bool reverse_matrix(void)
{
	int i, j, k, ik;
	memset(I, 0, sizeof(I));
	for(i = 1; i <= n; i++) I[i][i] = 1;
	for(k = 1; k <= n; k++)
	{
		double mkk, max = -inf;
		ik = -1;
		for(i = k; i <= n; i++)
		{
			if(a[i][k] > max)
			{
				max = a[i][k];
				ik = i;
			}
		}
		if(fabs(max) < eps) return false;
		if(ik != k)
		{
			for(i = 1; i <= n; i++)
			{
				swap(a[k][i], a[ik][i]);
				swap(I[k][i], I[ik][i]);
			}
		}
		for(i = 1; i <= n; i++)
		{
			if(i != k) a[i][k] /= -a[k][k];
			else mkk = 1 / a[k][k];
		}
		for(i = 1; i <= n; i++)
		{
			if(i != k)
			{
				for(j = k + 1; j <= n; j++) a[i][j] += a[i][k] * a[k][j];
				for(j = 1; j <= n; j++) I[i][j] += a[i][k] * I[k][j];
			}
		}
		for(j = k; j <= n; j++) a[k][j] *= mkk;
		for(j = 1; j <= n; j++) I[k][j] *= mkk;
	}
	return true;
}