矩陣求逆算法-全選主元高斯-約旦法

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: 全部代碼