C语言实现实数和复数矩阵及其各种运算（五）

一、前言

1. 本章开始，笔者将详细讲解矩阵的QR分解、特征值、特征向量等运算，并给出函数和测试demo的C代码，以及与matlab计算的结果；
2. 并且，本章相当于前面几章的大杂烩，前面所有的结构体、宏定义、函数本章基本全部都有用到。
3. 此外，由于本章求解QR矩阵/特征值/特征向量，实数和复数的过程原理完全一样，读者可以自行修改参数，简化复数计算的部分即可，或者完全不用修改，直接调用我的函数，虚部设置为0退化为实数运算即可；
4. 另，这将也是我此连载文章的系列精品之一，矩阵运算的功能越来越复杂，但是笔者全部分模块进行，化繁为简，化整为零，所以需要调用的小函数越来越多，读者遇到不熟悉的函数可以随时Jump到（一）、（二）、（三）、（四）的内容。
   （PS：另外这两个编辑器为啥功能不能优势互补，而且都没有贴matlab code的编辑器？）
   

二、QR分解

QR分解的数学原理读者可以直接点击链接复习一下，QR分解C语言实现我参考的是此大佬的。笔者这里整理出来了，C code直接给出：

/* QR Decompose */
void QR(const Matrix* A, Matrix* Q, Matrix* R)
{
   
	if (IsNullComplexMatrix(A))
	{
   
		puts("ERROE: An invalid matrix!\n");
		return;
	}
	// Not A Square Matrix
	if (A->row != A->column)
	{
   
		printf("ERROE: Not a square matrix!\n");
		return;
	}
	
	int i, j, k, m;
	int size;
	const int N = MatrixRow(A);
	ComplexType temp;

	// Column Vector Saving Column Vectors of A
	Matrix a, b;
	InitComplexMatrix(&a, N, 1);
	InitComplexMatrix(&b, N, 1);
	size = MatrixSize(A);
	if (MatrixSize(Q) != size)
	{
   
		// free(Q->arrayComplex);
		DestroyComplexMatrix(Q);                   // Free The Initial Matrix
		InitComplexMatrix(Q, A->row, A->column);   // Reset Size and Initialize
	}

	if (MatrixSize(R) != size)
	{
   
		// free(R->arrayComplex);
		DestroyComplexMatrix(R);
		InitComplexMatrix(R, A->row, A->column);
	}

	for (j = 0; j < N; ++j)
	{
   
		for (i = 0; i < N; ++i)
		{
   
			a.arrayComplex[i] = b.arrayComplex[i] = A->arrayComplex[i * A->column + j];   // Cols Vector of A
		}

		for (k = 0; k < j; ++k)
		{
   
			R->arrayComplex[k * R->column + j]._Val[0] = 0;
			R->arrayComplex[k * R->column + j]._Val[1] = 0;
//			InitComplex(R->arrayComplex + k * R->column + j);  // reset
			for (m = 0; m < N; ++m)
			{
   
				R->arrayComplex[k * R->column + j] = AddComplex(R->arrayComplex[k * R->column + j], \
					_Cmulcc(a.arrayComplex[m], Q->arrayComplex[m * Q->column + k]));
			}

			for (m = 0; m < N; ++m)
			{
   
				b.arrayComplex[m] = SubComplex(b.arrayComplex[m], _Cmulcc(R->arrayComplex[k * R->column + j], \
					Q->arrayComplex[m * Q->column + k]));
			}
		}

		temp = MatrixNorm2(&b);
		R->arrayComplex[j * R->column + j] = temp;

		for (i = 0; i < N; ++i)
		{
   
			Q->arrayComplex[i * Q->column + j] = DivComplex(b.arrayComplex[i], temp);
		}
	}
	// Free Local Memory
	Matrix ComplexMatrixArray[] = {
    a, b };
	int numDoubleMatrixArray = sizeof(ComplexMatrixArray) / sizeof(Matrix);
	DestroyComplexMatrixArray(ComplexMatrixArray, numDoubleMatrixArray);      // Complpex Matrix
    // OR
//	DestroyComplexMatrix(&a);
//	DestroyComplexMatrix(&b);
}

说明：
(1) 因为在后面的特征值的求解函数中需要用到，单单QR分解本身数学意义不大，因此单独的测试就免了；
(2) 后面释放结构体指针（也就是矩阵内存）的函数，调用的是DestroyComplexMatrixArray()而非DestroyComplexMatrix()，前面第一章给了二者实现的代码，读者留意一下，我个人建议读者定义的矩阵不多，最好不要调用这个destroy函数，老老实实一个一个调用DestroyComplexMatrix()/DestroyDoubleMatrix()函数释放内存；
(3) 另外，注意矩阵的个数numDoubleMatrixArray必须在调用DestroyComplexMatrixArray()函数之前求出来赋值给变量然后传给形参，而不能在DestroyComplexMatrixArray()函数里面求解，且ComplexMatrixArray必须定义为数组而不能是指针，因为无论是指针还是数组，实参传给形参时候，系统都会在函数内部将其当作指针处理，sizeof(ptr) ≡ 4Byte(64bit系统下)，而不是指向的那段内存的实际大小，但是sizeof(数组名)却是实际的那段内存的大小，这也是数组和指针的区别，C语言里面这一块也是一个比较容易出错的地方（笔者此处也差点搞忘了！！！）。

三、矩阵特征值

特征值计算比较复杂，关于特征值和特征向量的求解，线性代数理论知识大家可以参考此篇文章，关于特征值和特征向量C语言的实现，笔者是根据此位大佬修改的，这里我直接贴出笔者魔改后的代码：

/* eigen values */
void EigenValue(const Matrix* matrix, Matrix* eigenvalue)
{
   
	const int NUM = 100;   // Iteration Times
	
	// Local Matrice
	Matrix Q, R, temp;
	// Initiate
	InitComplexMatrix(&Q, matrix->row, matrix->column);
	InitComplexMatrix(&R, matrix->row, matrix->column);
	InitComplexMatrix(&temp, matrix->row, matrix->column);

	if (IsNullComplexMatrix(matrix))
	{
   
		puts("ERROE: An invalid matrix!\n");
		return;
	}
	// Copy matrix to temp in order not to change matrix
	CopyMatrix(matrix, &temp);
	// QR Decompose and 
	for (int k = 0; k < NUM; ++k)
	{
   
		QR(&temp, &Q, &R);
		MatrixMulMatrix(&R, &Q, &temp);
	}
	// Abstract Eigen Values from the Diagonal Elements of temp =  Q * R
	for (int k = 0; k < temp.row; ++k)
	{
   
		eigenvalue->arrayComplex[k] = temp.arrayComplex[k * temp.column + k];
	}
	// Free Local Memory
	Matrix ComplexMatrixArray[] = {
    R, Q, temp };
	int numComplexMatrixArray = sizeof(ComplexMatrixArray) / sizeof(Matrix);
	DestroyComplexMatrixArray(ComplexMatrixArray, numComplexMatrixArray);      // Complpex Matrix
//	return eigenvalue;
}

说明：
(1) 因为在后面的特征向量的求解中需要用到，因此单独的测试就免了；
(2) 同样这里销毁结构体内部指针（即矩阵）内存调用的也是一次销毁函数DestroyComplexMatrixArray()；
(3) 读者注意一下笔者在此特征值求解函数中调用的一些函数，如果忘记了可以跳到前面几章节review。

三、矩阵特征向量

我是一个没有感情的复读机：特征向量的计算更为复杂，关于特征向量的求解，线性代数理论知识大家可以参考此篇文章，关于特征向量C语言的实现，笔者是根据此位大佬修改的，这里我直接贴出笔者修魔改后的代码：

/* Negative Complex : Complex_B = -Complex_A = -creal(Complex_A) - cimag(Complpex_A) */
ComplexType NegativeComplex(const ComplexType Complex_A)
{
   
	ComplexType Complex_B;
	Complex_B = _Cmulcr(Complex_A, -1.0);
	// OR
/*
	Complex_B._Val[0] = -creal(Complex_A);
	Complex_B._Val[1] = cimag(Complex_A) * (-1.0);
*/
	return Complex_B;
}

/* eigen vectors */
void EigenVector(const Matrix* matrix, const Matrix* eigenvalue, Matrix* eigenvector)
{
   
	if (IsNullComplexMatrix(matrix) || IsNullComplexMatrix(eigenvalue))
	{
   
		puts("ERROE: An invalid matrix!\n");
		return;
	}
	int i, j, q;
	int m;

	// Access to Eigen Values
	int count;
	int num = MatrixRow(matrix);   // = matrix->row: Numbers of Eigen Values or Cols
	ComplexType evalue;

	// Access to temp
	ComplexType sum, midsum, mid;
	Matrix temp;   // temp = A - λI: (A - λI) * x = 0
	InitComplexMatrix(&temp, matrix->row, matrix->column);

	for (count = 0; count < num; ++count)
	{
   
		// Calculate x: Ax = λ * x
		evalue = eigenvalue->arrayComplex[</