C语言实现复矩阵求逆（四阶）

1. 伴随矩阵和代数余子式法
2. 高斯消元法
3. LU分解

前言： 接上一篇：实现矩阵的求逆
参考：
代码实现矩阵求逆的三种方式
矩阵的行列式
复数矩阵计算行列式

1. 伴随矩阵和代数余子式法

1.1 原理

见上篇

1.2 代码

#include <stdio.h>
#include <math.h>
typedef unsigned char uint8;
typedef struct {
	double re;
	double im;
} Complex;

//全局变量定义
uint8 MatrixSize = 4;
Complex MatrixTest[4][4] = { 0 };
Complex AdjMatrix[4][4] = { 0 };
Complex detA4 = { 0 };
Complex InvMtrx[4][4] = { 0 };


Complex GetdetA3(Complex matrix[3][3]);
void Adj_matrix4(Complex* matrix);
void inv(Complex* matrixtest);
Complex add(Complex a, Complex b);
Complex sub(Complex a, Complex b);
Complex Mul(Complex a, Complex b);
Complex Mul3(Complex a, Complex b, Complex c);
Complex Div(Complex a, Complex b);


int  main()
{
	uint8 i = 0,j=0;
	printf("Hello!");
	//matlab随机生成满秩矩阵
	Complex matrixtest[4][4] = { {{0.7094,0.8909},{0.7547,0.9593},{0.2760,0.5472},{0.6797,0.1386}},
		{{0.6551,0.1493},{0.1626,0.2575},{0.1190,0.8407},{0.4984,0.2543}},
		{{0.9597,0.8143},{0.3404,0.2435},{0.5853,0.9293},{0.2238,0.3500}},
		{{0.7513,0.1966},{0.2551,0.2511},{0.5060,0.6160},{0.6991,0.4733}} };
	//double det = 0;
	inv(matrixtest);
	return 0;
}

void inv(Complex* matrixtest) {

	if (matrixtest == NULL)   /* 检查输入的指针是否为空*/
	{
		printf("matrix pointer is NULL.\n");
	}
	Adj_matrix4(matrixtest);

	if (detA4.re != 0 && detA4.im != 0) {
		for (uint8 i = 0; i < 4; i++) {
			for (uint8 j = 0; j < 4; j++) {
				InvMtrx[i][j] = Div(AdjMatrix[i][j], detA4);
			}
		}
	}
}

Complex GetdetA3(Complex matrix[3][3])
{
	Complex detA3 = { 0 };
	Complex a1 = { 0 }, a2 = { 0 }, a3 = { 0 };
	Complex b1 = { 0 }, b2 = { 0 }, b3 = { 0 };
	a1 = Mul3(matrix[0][0], matrix[1][1], matrix[2][2]);
	a2 = Mul3(matrix[0][1], matrix[1][2], matrix[2][0]);
	a3 = Mul3(matrix[0][2], matrix[1][0], matrix[2][1]);
	b1 = Mul3(matrix[0][2], matrix[1][1], matrix[2][0]);
	b2 = Mul3(matrix[0][1], matrix[1][0], matrix[2][2]);
	b3 = Mul3(matrix[0][0], matrix[1][2], matrix[2][1]);
	detA3 = sub(add(add(a1, a2), a3), add(add(b1, b2), b3));
	return detA3;
}

//四阶矩阵的伴随矩阵
void Adj_matrix4(Complex* matrix) {
	Complex tmp[3][3] = { {0} ,{0} ,{0} };
	for (uint8 row = 0; row < 4; row++) {
		for (uint8 colm = 0; colm < 4; colm++) {
			//求第row、colm余子式矩阵
			uint8 k = 0, p = 0;
			for (uint8 i = 0; i < 3; i++) {
				if (i == colm) p++;
				for (uint8 j = 0; j < 3; j++) {
					if (j == row) k++;
					tmp[i][j] = *(matrix + k + 4 * p + j);
				}
				k = 0;
				p++;
			}
			//求余子式矩阵行列式
			Complex C_0colm = GetdetA3(tmp);
			//求出相应的代数余子式
			AdjMatrix[row][colm].re = pow(-1, row + colm) * C_0colm.re;
			AdjMatrix[row][colm].im = pow(-1, row + colm) * C_0colm.im;
			if (row == 0) detA4 = add(detA4, Mul((*(matrix + row + colm * 4)), AdjMatrix[row][colm]));
		}
	}
}



Complex add(Complex a, Complex b)
{
	Complex c = { 0 };
	c.re = a.re + b.re;
	c.im = a.im + b.im;
	return c;
}
Complex sub(Complex a, Complex b)
{
	Complex c = { 0 };
	c.re = a.re - b.re;
	c.im = a.im - b.im;
	return c;
}
Complex Mul(Complex a, Complex b)
{
	Complex c = { 0 };
	c.re = a.re * b.re - a.im * b.im;
	c.im = a.re * b.im + b.re * a.im;
	return c;
}
Complex Mul3(Complex a, Complex b, Complex c)
{
	Complex d = { 0 };
	d = Mul(a, b);
	d = Mul(d, c);
	return d;
}
Complex Div(Complex a, Complex b)
{
	Complex c = { 0 };
	c.re = (a.re * b.re + a.im * b.im) / (b.re * b.re + b.im * b.im);
	c.im = (a.im * b.re - a.re * b.im) / (b.re * b.re + b.im * b.im);
	return c;
}