矩阵求逆算法-全选主元高斯-约旦法

矩阵求逆算法-全选主元高斯-约旦法

Tags： 逆矩阵

---

全选主元高斯-约旦法求逆的步骤如下：

1. 对于 k 从 0 到 n - 1 作如下几步：

1. 从第 k 行、第 k 列开始的右下角子阵中选取绝对值最大的元素，并记住次元素所在的行号和列号，在通过行交换和列交换将它交换到主元素位置上。这一步称为全选主元。
2. m(k, k) = 1 / m(k, k)
3. m(k, j) = m(k, j) * m(k, k)，j = 0, 1, ..., n-1；j != k
4. m(i, j) = m(i, j) - m(i, k) * m(k, j)，i, j = 0, 1, ..., n-1；i, j != k
5. m(i, k) = -m(i, k) * m(k, k)，i = 0, 1, ..., n-1；i != k

2. 根据在全选主元过程中所记录的行、列交换的信息进行恢复，恢复的原则如下：

在全选主元过程中，先交换的行（列）后进行恢复；原来的行（列）交换用列（行）交换来恢复。

---

示例代码

1. 矩阵结构定义

#ifdef DEBUG
#define LM_ASSERT(x)	do{ if(!(x)) while(1); }while(0)
#else
#define LM_ASSERT(x)
#endif

typedef float real_t;

typedef struct matrix
{
	unsigned int row;
	unsigned int col;
	real_t *m;
}matrix_t;

#define MATRIX(M, r, c) (*((M)->m + (r)*(M)->col + (c)))

2 函数实现

#define FABS(x)		fabs(x)

/**
 * @brief 交换行
 */
void matrix_swap_row(matrix_t *m, unsigned int i, unsigned int j)
{
	unsigned int k;
	real_t tmp;

	LM_ASSERT(i < m->row);
	LM_ASSERT(j < m->row);

	for(k=0; k<m->col; k++)
	{
		tmp = MATRIX(m, i, k);
		MATRIX(m, i, k) = MATRIX(m, j, k);
		MATRIX(m, j, k) = tmp;
	}
}

/**
 * @brief 交换列
 */
void matrix_swap_col(matrix_t *m, unsigned int i, unsigned int j)
{
	unsigned int k;
	real_t tmp;

	LM_ASSERT(i < m->col);
	LM_ASSERT(j < m->col);

	for(k=0; k<m->row; k++)
	{
		tmp = MATRIX(m, k, i);
		MATRIX(m, k, i) = MATRIX(m, k, j);
		MATRIX(m, k, j) = tmp;
	}
}

/**
 * @brief 复制矩阵
 */
void matrix_copy(matrix_t *to, matrix_t *from)
{
	unsigned int i, j;

	matrix_reshape(to, from->row, from->col);

	for(i=0; i<from->row; i++)
		for(j=0; j<from->col; j++)
			MATRIX(to, i, j) = MATRIX(from, i, j);
}

/**
 * @brief 设置矩阵大小
 */
int matrix_reshape(matrix_t *m, unsigned int row, unsigned int col)
{
	LM_ASSERT(m != NULL);
	LM_ASSERT(row != 0);
	LM_ASSERT(col != 0);

	if (row * col == m->row * m->col)
	{
		m->row = row;
		m->col = col;
	}
	else
	{
		if (m->m != NULL)
			free(m->m);

		m->m = malloc(row * col * sizeof(real_t));
		if (m->m != NULL)
		{
			m->row = row;
			m->col = col;
		}
		else
		{
			m->row = m->col = 0;
			return -1;
		}
	}

	return 0;
}

/**
 * @brief 求逆矩阵
 */
int matrix_inv(matrix_t *inv, matrix_t *a)
{
	int i, j, k;
	int ret = 0;

	//! 必须是方阵
	LM_ASSERT(a->row == a->col);

	unsigned int is[a->row];
	unsigned int js[a->col];

	real_t max;

	matrix_reshape(inv, a->row, a->col);

	matrix_copy(inv, a);

	for(k=0; k<inv->row; k++)
	{
		//step 1, 全选主元
		max = 0;
		is[k] = k;
		js[k] = k;

		for(i=k; i<inv->row; i++)
		{
			for(j=k; j<inv->col; j++)
			{
				if(max < FABS(MATRIX(inv, i, j)))
				{
					max = FABS(MATRIX(inv, i, j));
					is[k] = i;
					js[k] = j;
				}
			}
		}

		if(max < 0.0001)
		{	//! 无逆矩阵
			ret = -1;
			goto end;
		}

		//交换
		if(is[k] != k)
		{
			matrix_swap_row(inv, k, is[k]);
		}
		if(js[k] != k)
		{
			matrix_swap_col(inv, k, js[k]);
		}

		MATRIX(inv, k, k) = 1 / MATRIX(inv, k, k);

		for(j=0; j<inv->col; j++)
		{
			if(j != k)
				MATRIX(inv, k, j) *= MATRIX(inv, k, k);
		}
		for(i=0; i<inv->row; i++)
		{
			if(i != k)
			{
				for(j=0; j<inv->col; j++)
				{
					if(j != k)
						MATRIX(inv, i, j) -= MATRIX(inv, i, k) * MATRIX(inv, k, j);
				}
			}
		}
		for(i=0; i<inv->row; i++)
		{
			if(i != k)
				MATRIX(inv, i, k) *= -MATRIX(inv, k, k);
		}

	}

	//恢复
	//本来 row <-> is[k], column <-> js[k]
	//恢复时：row <-> js[k], column <-> is[k]
	for(k=inv->row-1; k>=0; k--)
	{
		if(js[k] != k)
		{
			matrix_swap_row(inv, k, js[k]);
		}
		if(is[k] != k)
		{
			matrix_swap_col(inv, k, is[k]);
		}
	}

end:
	return ret;;
}

---

GitHub: 全部代码