一、前言

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[count];
		CopyMatrix(matrix, &temp);
		for (i = 0; i < temp.column; ++i)
		{
			temp.arrayComplex[i * temp.column + i] = SubComplex(temp.arrayComplex[i * temp.column + i], evalue);
			//			temp->arrayComplex[i * temp->column + i]._Val[0] -= creal(evalue);
			//			temp->arrayComplex[i * temp->column + i]._Val[1] -= cimag(evalue);
		}

		// Transform temp to Ladder Matrix
		for (i = 0; i < temp.row - 1; ++i)
		{
			mid._Val[0] = creal(temp.arrayComplex[i * temp.column + i]);   // Diagonal Element
			mid._Val[1] = cimag(temp.arrayComplex[i * temp.column + i]);
			for (j = i; j < temp.column; ++j)
			{
				temp.arrayComplex[i * temp.column + j] = DivComplex(temp.arrayComplex[i * temp.column + j], mid);
			}

			for (j = i + 1; j < temp.row; ++j)
			{
				mid = temp.arrayComplex[j * temp.column + i];
				for (q = i; q < temp.column; ++q)
				{
					temp.arrayComplex[j * temp.column + q] = SubComplex(temp.arrayComplex[j * temp.column + q], \
						_Cmulcc(temp.arrayComplex[i * temp.column + q], mid));

				}
			}
		}
		midsum._Val[0] = 1; 
		midsum._Val[1] = 0;
		eigenvector->arrayComplex[(eigenvector->row - 1) * eigenvector->column + count]._Val[0] = 1;
		eigenvector->arrayComplex[(eigenvector->row - 1) * eigenvector->column + count]._Val[1] = 0;
		for (m = temp.row - 2; m >= 0; --m)
		{
//			InitComplex(&sum);
			sum._Val[0] = 0; sum._Val[1] = 0;   // Zero Complex
			for (j = m + 1; j < temp.column; ++j)
			{
				sum = AddComplex(sum, _Cmulcc(temp.arrayComplex[m * temp.column + j],
					eigenvector->arrayComplex[j * eigenvector->column + count]));
			}
			sum = DivComplex(NegativeComplex(sum), temp.arrayComplex[m * temp.column + m]);  // Warning: Parameters' Type
			//sum = -sum / *(temp.arrayComplex[m * temp.column + m]);
			midsum = AddComplex(midsum, _Cmulcc(sum, sum));
			eigenvector->arrayComplex[m * eigenvector->column + count] = sum;
		}

		midsum = csqrt(midsum);
		for (i = 0; i < eigenvector->row; ++i)  // One Column Vector--Eigen Vector
		{
			eigenvector->arrayComplex[i * eigenvector->column + count] = 
				DivComplex(eigenvector->arrayComplex[i * eigenvector->column + count], midsum);
		}
	}
	DestroyComplexMatrix(&temp);
//	return eigenvector;
}

说明：
(1) 第一个函数用于求解复数的相反数，即a+bi—> -a-bi，当然也可以直接调用complex.h头文件里的_Cmulcr()函数：

sum = _Cmulcr(sum, -1.0);

(2) 这里只有一个local matrix variable，因此笔者直接调用的DestroyComplexMatrix()函数销毁内存；
(3) 笔者注释一定要看，包括前面的特征值计算，这样对比数学原理，函数理解起来非常简单。
(4) 将本章的QR分解、特征值、特征向量函数放在一起测试，测试的demo如下：

#include<stdio.h>
#include<stdlib.h>
#include<complex.h>
#include <math.h>

typedef _Dcomplex ComplexType;
typedef double DoubleType;

typedef struct
{
	int row, column;
	ComplexType* arrayComplex;
}Matrix;

typedef struct
{
	int row, column;
	DoubleType* arrayDouble;
}Matrix2Double;

typedef enum
{
	False = 0, True = 1
}Bool;

/*
  Initiation of Complex Number
 */
void InitComplex(ComplexType* Complex)
{
	Complex->_Val[0] = 0.0;
	Complex->_Val[1] = 0.0;

}

/*
  Validity of Complex Matrix
*/
Bool IsNullComplexMatrix(const Matrix* matrix)
{
	int size = matrix->row * matrix->column;

	if (size <= 0 || matrix->arrayComplex == NULL)
	{
		return True;
	}
	return False;
}


/*
  Initiation of Complex Matrix
 */
void InitComplexMatrix(Matrix* matrix, int row, int column)
{
	int size = row * column * sizeof(ComplexType);
	if (size <= 0)
	{
		puts("ERROE: An invalid matrix!\n");
		return;
	}
	matrix->arrayComplex = (ComplexType*)malloc(size);
	if (matrix->arrayComplex)
	{
		matrix->row = row;
		matrix->column = column;
		for (int row_i = 0; row_i < row; ++row_i)
		{
			for (int column_j = 0; column_j < column; ++column_j)
			{
				InitComplex(matrix->arrayComplex + row_i * matrix->column + column_j);

			}
		}

	}
}

/*
  Free Memory of Complex Matrix
*/
void DestroyComplexMatrix(Matrix* matrix)
{
	if (!IsNullComplexMatrix(matrix))
	{
		free(matrix->arrayComplex);
		matrix->arrayComplex = NULL;
	}
	matrix->row = matrix->column = 0;
}

/*
  Free Memory of Complex Matrice Array
*/
void DestroyComplexMatrixArray(Matrix matrixArray[], int num)  // Array Transfer--->Pointer Transfer
{
	if (num)      // if no cell
	{
		for (int i = 0; i < num; i++)
		{
			DestroyComplexMatrix(&matrixArray[i]);  // Nested Call of DestroyComplexMatrix()
		}
	}
}

/*
  return matrix row size
*/
int MatrixRow(const Matrix* matrix)
{
	return matrix->row;
}

/*
  return matrix column size
*/
int MatrixColumn(const Matrix* matrix)
{
	return matrix->column;
}

/*
  return matrix column size
*/
int MatrixSize(const Matrix* matrix)
{
	return (matrix->row * matrix->column);
}

/*
   Add Complex: Complex_C = Complex_A + Complex_B
*/
ComplexType AddComplex(const ComplexType Complex_A, const ComplexType Complex_B)
{
	ComplexType Complex_C;
	Complex_C._Val[0] = creal(Complex_A) + creal(Complex_B);
	Complex_C._Val[1] = cimag(Complex_A) + cimag(Complex_B);
	return Complex_C;
}

/*
   Subvision Complex: Complex_C = Complex_A - Complex_B
*/
ComplexType SubComplex(const ComplexType Complex_A, const ComplexType Complex_B)
{
	ComplexType Complex_C;
	Complex_C._Val[0] = creal(Complex_A) - creal(Complex_B);
	Complex_C._Val[1] = cimag(Complex_A) - cimag(Complex_B);
	return Complex_C;
}

/*
   Division Complex : Complex_C = Complex_A / Complex_B
*/
ComplexType DivComplex(const ComplexType Complex_A, const ComplexType Complex_B)   // Tip: Return Value Is NOT A Pointer
{
	ComplexType Complex_C;
	Complex_C._Val[0] = (creal(Complex_A) * creal(Complex_B) + cimag(Complex_A)
		* cimag(Complex_B)) / (pow(creal(Complex_B), 2) + pow(cimag(Complex_B), 2));
	Complex_C._Val[1] = (cimag(Complex_A) * creal(Complex_B) - creal(Complex_A)
		* cimag(Complex_B)) / (pow(creal(Complex_B), 2) + pow(cimag(Complex_B), 2));
	return Complex_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);
	return Complex_B;
}

/*
   2-norm of a Matrix
*/
ComplexType MatrixNorm2(const Matrix* matrix)
{
	ComplexType norm;
	norm._Val[0] = 0; norm._Val[1] = 0;
	if (IsNullComplexMatrix(matrix))
	{
		puts("ERROE: An invalid matrix!\n");
		return (csqrt(norm));
	}
	else
	{
		for (int i = 0; i < matrix->row; ++i)
		{
			for (int j = 0; j < matrix->column; ++j)
				norm = AddComplex(norm, _Cmulcc(matrix->arrayComplex[i * matrix->column + j], matrix->arrayComplex[i * matrix->column + j]));
		}
		return (csqrt(norm));
	}
}

/*
  Copy: Complex matrixA = Complex matrixB
*/
void CopyMatrix(const Matrix* matrix_A, Matrix* matrix_B)
{
	if (IsNullComplexMatrix(matrix_A))
	{
		puts("ERROE: An invalid matrix!\n");
		return;
	}
	else
	{
		int index = 0;
		for (int row_i = 0; row_i < matrix_A->row; row_i++)
		{
			for (int column_j = 0; column_j < matrix_A->column; column_j++)
			{
				index = matrix_B->column * row_i + column_j;
				matrix_B->arrayComplex[index] = matrix_A->arrayComplex[index];
			}
		}
	}
}

/*
   QR Decompose
*/
void QR(const Matrix* A, Matrix* Q, Matrix* R)
{
	if (IsNullComplexMatrix(A))
	{
		puts("ERROE: An invalid matrix!\n");
		return;
	}
	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;

	Matrix a, b;
	InitComplexMatrix(&a, N, 1);
	InitComplexMatrix(&b, N, 1);
	size = MatrixSize(A);
	if (MatrixSize(Q) != size)
	{
		DestroyComplexMatrix(Q);                   
		InitComplexMatrix(Q, A->row, A->column); 
	}

	if (MatrixSize(R) != size)
	{
		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];
		}

		for (k = 0; k < j; ++k)
		{
			R->arrayComplex[k * R->column + j]._Val[0] = 0;
			R->arrayComplex[k * R->column + j]._Val[1] = 0;
			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);
		}
	}
	Matrix ComplexMatrixArray[] = { a, b };
	int numDoubleMatrixArray = sizeof(ComplexMatrixArray) / sizeof(Matrix);
	DestroyComplexMatrixArray(ComplexMatrixArray, numDoubleMatrixArray);
}

/*
  Complex Matrix Multiple: matrixC = matrixA * matrixB
*/
void MatrixMulMatrix(const Matrix* matrixA, const Matrix* matrixB, Matrix* matrixC)
{
	if (IsNullComplexMatrix(matrixA) || IsNullComplexMatrix(matrixB))
	{
		puts("ERROE: An invalid matrix!\n");
		return;
	}

	else if (matrixA->column != matrixB->row)
	{
		puts("ERROE: An incompatable matrix!\n");
		return;
	}
	else
	{
		int row_i, column_j, ij;
		int indexA, indexB, indexC;
		ComplexType tempCell_C;

		for (row_i = 0; row_i < matrixC->row; ++row_i)
		{
			for (column_j = 0; column_j < matrixC->column; ++column_j)
			{
				tempCell_C._Val[0] = 0;
				tempCell_C._Val[1] = 0;
				for (ij = 0; ij < matrixA->column; ++ij)
				{
					indexA = row_i * matrixA->column + ij;
					indexB = ij * matrixB->column + column_j;
					tempCell_C = AddComplex(tempCell_C,
						_Cmulcc(matrixA->arrayComplex[indexA], matrixB->arrayComplex[indexB]));
				}
				indexC = row_i * matrixC->column + column_j;
				matrixC->arrayComplex[indexC] = tempCell_C;
			}
		}
	}
}

/*
   eigen values
*/

void EigenValue(const Matrix* matrix, Matrix* eigenvalue)
{
	const int NUM = 100; 

	Matrix Q, R, temp;
	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;
	}

	CopyMatrix(matrix, &temp);

	for (int k = 0; k < NUM; ++k)
	{
		QR(&temp, &Q, &R);
		MatrixMulMatrix(&R, &Q, &temp);
	}
	for (int k = 0; k < temp.row; ++k)
	{
		eigenvalue->arrayComplex[k] = temp.arrayComplex[k * temp.column + k];
	}
	Matrix ComplexMatrixArray[] = { R, Q, temp };
	int numComplexMatrixArray = sizeof(ComplexMatrixArray) / sizeof(Matrix);
	DestroyComplexMatrixArray(ComplexMatrixArray, numComplexMatrixArray);
}

/*
   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;
	int count;
	int num = MatrixRow(matrix); 
	ComplexType evalue;

	ComplexType sum, midsum, mid;
	Matrix temp;
	InitComplexMatrix(&temp, matrix->row, matrix->column);

	for (count = 0; count < num; ++count)
	{
		evalue = eigenvalue->arrayComplex[count];
		CopyMatrix(matrix, &temp);
		for (i = 0; i < temp.column; ++i)
		{
			temp.arrayComplex[i * temp.column + i] = SubComplex(temp.arrayComplex[i * temp.column + i], evalue);
		}
		for (i = 0; i < temp.row - 1; ++i)
		{
			mid._Val[0] = creal(temp.arrayComplex[i * temp.column + i]); 
			mid._Val[1] = cimag(temp.arrayComplex[i * temp.column + i]);
			for (j = i; j < temp.column; ++j)
			{
				temp.arrayComplex[i * temp.column + j] = DivComplex(temp.arrayComplex[i * temp.column + j], mid);
			}

			for (j = i + 1; j < temp.row; ++j)
			{
				mid = temp.arrayComplex[j * temp.column + i];
				for (q = i; q < temp.column; ++q)
				{
					temp.arrayComplex[j * temp.column + q] = SubComplex(temp.arrayComplex[j * temp.column + q], \
						_Cmulcc(temp.arrayComplex[i * temp.column + q], mid));

				}
			}
		}
		midsum._Val[0] = 1;
		midsum._Val[1] = 0;
		eigenvector->arrayComplex[(eigenvector->row - 1) * eigenvector->column + count]._Val[0] = 1;
		eigenvector->arrayComplex[(eigenvector->row - 1) * eigenvector->column + count]._Val[1] = 0;
		for (m = temp.row - 2; m >= 0; --m)
		{
			sum._Val[0] = 0; sum._Val[1] = 0;
			for (j = m + 1; j < temp.column; ++j)
			{
				sum = AddComplex(sum, _Cmulcc(temp.arrayComplex[m * temp.column + j],
					eigenvector->arrayComplex[j * eigenvector->column + count]));
			}
			sum = DivComplex(NegativeComplex(sum), temp.arrayComplex[m * temp.column + m]);
			midsum = AddComplex(midsum, _Cmulcc(sum, sum));
			eigenvector->arrayComplex[m * eigenvector->column + count] = sum;
		}

		midsum = csqrt(midsum);
		for (i = 0; i < eigenvector->row; ++i)
		{
			eigenvector->arrayComplex[i * eigenvector->column + count] =
				DivComplex(eigenvector->arrayComplex[i * eigenvector->column + count], midsum);
		}
	}
	DestroyComplexMatrix(&temp);
}

int main(void)
{
	Matrix matrix;
	InitComplexMatrix(&matrix, 3, 3);  // 3 * 3	
	matrix.arrayComplex[0]._Val[0] = 2;  matrix.arrayComplex[0]._Val[1] = 2;
	matrix.arrayComplex[1]._Val[0] = 3;  matrix.arrayComplex[1]._Val[1] = 7;
	matrix.arrayComplex[2]._Val[0] = 3;  matrix.arrayComplex[2]._Val[1] = -6;
	matrix.arrayComplex[3]._Val[0] = 6;  matrix.arrayComplex[3]._Val[1] = 2;
	matrix.arrayComplex[4]._Val[0] = -7; matrix.arrayComplex[4]._Val[1] = 7;
	matrix.arrayComplex[5]._Val[0] = 2;  matrix.arrayComplex[5]._Val[1] = 4;
	matrix.arrayComplex[6]._Val[0] = -9;  matrix.arrayComplex[6]._Val[1] = 6;
	matrix.arrayComplex[7]._Val[0] = 5;  matrix.arrayComplex[7]._Val[1] = -4;
	matrix.arrayComplex[8]._Val[0] = -2;  matrix.arrayComplex[8]._Val[1] = -1;

	Matrix EigenValues;
	InitComplexMatrix(&EigenValues, MatrixRow(&matrix), 1);  // 3 * 1
	EigenValue(&matrix, &EigenValues);

	Matrix EigenVectors;
	InitComplexMatrix(&EigenVectors, MatrixRow(&matrix), MatrixColumn(&matrix));  // 3 * 3
	EigenVector(&matrix, &EigenValues, &EigenVectors);

	for (int i = 0; i < EigenValues.row; i++)
		for (int j = 0; j < EigenValues.column; j++)
			printf("%lf, %lfi\n", creal(EigenValues.arrayComplex[i * EigenValues.column + j]), cimag(EigenValues.arrayComplex[i * EigenValues.column + j]));

	printf("\n\n\n\n\n");

	for (int index = 0; index < MatrixSize(&EigenVectors); index++)
		printf("%lf, %lfi\n", creal(EigenVectors.arrayComplex[index]), cimag(EigenVectors.arrayComplex[index]));
	// Free Array of Memories of Complex Matrice by Calling ComplexMatrixArray() 
	Matrix ComplexMatrixArray[] = { matrix, EigenValues, EigenVectors };
	int numComplexMatrixArray = sizeof(ComplexMatrixArray) / sizeof(Matrix);
	DestroyComplexMatrixArray(ComplexMatrixArray, numComplexMatrixArray);
	
	return 0;

测试结果我直接通过截图的形式给出，VS C（上）和matlab（下）：


PS：
(1) <2020.8.20补充> 上面看起来，二者的特征向量不一样，我刚开始也有点纳闷，后来经小伙伴@hamsterWalle的提示，发现matlab调用eig()函数求导的特征向量并未对向量进行归一化（归一化：特征向量/列向量各个元素的平方和为1），而在笔者的C程序中，对其进行了归一化，因此结果不一样，这里我做了两组实验：
实验一：注意到前面给的求特征向量的函数EigenVector()中有两句程序

midsum = AddComplex(midsum, _Cmulcc(sum, sum));

midsum = csqrt(midsum);   // complex.h: csqrt()

即是对特征向量进行归一化处理。此外，为了验证此函数求得的特征向量确实是归一化之后的，笔者在上面的demo中，main()函数之前加了一个归一化函数（注意是列向量，因此列向量的元素index与行向量有区别）：

/*
	Normalization of Column Vector
*/
void NormVector(const Matrix* EigenVector, Matrix* NormEigenVector)
{
	if (IsNullComplexMatrix(EigenVector))
		return;
	else if ((MatrixRow(EigenVector) != MatrixRow(NormEigenVector)) || (MatrixColumn(EigenVector) != MatrixColumn(NormEigenVector)))
		return;
	else
	{
		ComplexType exp = { 2, 0 };  // initiate a complex number
		for (int column_j = 0; column_j < EigenVector->column; ++column_j)
		{
			ComplexType magnitude; InitComplex(&magnitude);
			for (int row_i = 0; row_i < EigenVector->row; row_i++)
			{
				int index = row_i * EigenVector->column + column_j;
				magnitude = AddComplex(magnitude, cpow(EigenVector->arrayComplex[index], exp));
// OR: cpow() -----> _Cmulcc()
//				magnitude = AddComplex(magnitude,  _Cmulcc(EigenVector->arrayComplex[index], EigenVector->arrayComplex[index]));
			}

			magnitude = csqrt(magnitude);

			for (int row_i = 0; row_i < EigenVector->row; row_i++)
			{
				int index = row_i * NormEigenVector->column + column_j;
				NormEigenVector->arrayComplex[index] = DivComplex(EigenVector->arrayComplex[index], magnitude);
			}
		}	
	}
}

另外main函数内加上：

	Matrix NormEigenVectors;
	InitComplexMatrix(&NormEigenVectors, MatrixRow(&EigenVectors), MatrixColumn(&EigenVectors));  // 3 * 3
	NormVector(&EigenVectors, &NormEigenVectors);
// AND
		printf("\n\n\n\n\n");
	for (int index = 0; index < MatrixSize(&NormEigenVectors); index++)
		printf("%lf, %lfi\n", creal(NormEigenVectors.arrayComplex[index]), cimag(NormEigenVectors.arrayComplex[index]));

此外，如果要释放NormEIgenVectors矩阵的内存（其实是结构体内部指针的内存），还需要加一笔：

Matrix ComplexMatrixArray[] = { matrix, EigenValues, EigenVectors, NormEigenVectors };

最后编译结果完全一样，VS C结果如下，说明EigenVector()求得的特征向量已经归一化了，不需要单独写一个函数归一化：

实验二：为了验证上面VS和matlab求得的特征向量结果本质完全一样，只是matlab的eig()函数没有对特征向量归一化，笔者这里用matlab手动对上面matlab计算结果进行归一化，以特征值-12.060159-4.426611i对应的特征向量为例：

matlab结果截图如下：


与前面VS求得的第一个特征值对应的特征向量（如下）完全一样：

(2) 关于matlab的eig()函数，若返回两个矢量，第一个是特征向量组成的矩阵，第二个是特征值组成的对角矩阵；若默认返回一个参数，则是一个特征值组成的列向量，此时matlab计算结果如下：

(3) 在求解特征向量函数中调用了前面章节没有涉及到的新的complex.h库函数—csqrt()函数和后面的归一化函数NormVector()中的cpow()函数，函数原型分别如下：

_ACRTIMP _Dcomplex __cdecl csqrt(_In_ _Dcomplex _Z);

_ACRTIMP _Dcomplex __cdecl cpow(_In_ _Dcomplex _X, _In_ _Dcomplex _Y);

(4) 无论是matlab调用的eig()函数还是笔者编写的EigenValue()和EigenVector()函数，请读者注意维度，另外，求得的特征值和特征向量（矩阵）都是以列向量依次对应书写的，笔者都有标出。
(5) <2020.8.20补充> 特征值若是实数，调用matlab的eig()函数会自动给特征值从小到大进行排序，读者若想实现跟matlab一样的效果，可以再写一个排序函数（快速/冒泡等排序算法），在调用完EigenValue()求得乱序的特征值之后，再调用该排序函数对特征值从小到大排序，然后再调用EigenValue()函数求特征向量，这样求得的特征值和特征向量会与matlab求得的顺序完全一致。

四、总结

实数/复数矩阵的常见运算就到此为止，前后一起花了我大概四天的时间完成了这五篇连载博客的撰写，写demo、测试、总结等，希望能帮助到各位社区小伙伴！

如何理解特征值为复数的情况

特征值与特征向量

特征值可定义为，若有$Ax=\lambda x$，则称$x$为$A$的特征向量，$\lambda$为相应的特征值。这时我们可以发现，如果$\lambda$是实数，那么矩阵$A$对向量$x$的作用，恰恰是将向量$x$放大或缩小了一定倍数，而方向不变。如果特征值绝对值大于1，那么$x$的长度将被放大（方向可能改变），如果特征值绝对值小于1，那么$x$的长度被缩小（方向也可能改变）。

特征向量在控制领域的意义

由于答主本科阶段学的是自动化专业，接下来用控制领域的知识进一步阐述特征向量相关知识。特征向量是什么呢？就是系统状态的某个特殊值。而根据线性定理，总可以将系统状态拆成不同特征向量的线性组合，所以也可以认为特征向量是系统状态的一个分量。在（离散）控制系统中，若$\varphi$是系统离散域的状态转移矩阵，假如某个时刻系统的状态$x$，恰好等于特征向量的倍数，那么每过一拍后，新的状态$x(n+1)=\varphi x(n)=\lambda x(n)$。如果$\lambda$的绝对值大于1，那么$x$中的值将随时间推移而不断变大，如果$\lambda$绝对值小于1，那么$x$中的值将随时间推移而不断衰减。系统的输出也包括在系统状态中，因此如果$\lambda$绝对值大于1，系统就是发散的。这也正是Z变换中的结论。

特征值为复数的情况

题主问的是特征值为复数的情况。当$\lambda$为复数时，从数学的角度看，可以认为状态$x$每经过一个控制节拍（经过一拍意味着系统经过了一次状态转移，下同），不仅长度发生改变，还产生了一个复平面内的旋转。那么复平面上的旋转，与实际的物理意义有什么对应关系呢？在我看来，可以这么理解：若选取合适坐标系，复平面中的值，其实部表征了当前物理量的大小，虚部象征了物理量转化（如能量转移）的趋势。以理想的LC串联电路为例：

电路工作过程为：若电感的感应电动势大于电容，则它给电容充能，否则电容给电感充能。由于没有能量的损耗，这个过程会持续到天荒地老。因此，若选取电容电压和电感电流构成一组系统状态，其离散系统状态转移矩阵的特征值就是“绝对值为1，而角度不为0”的一个复数。

举个例子，若电容为$1F$，电感为$1H$，系统的初始状态（电容电压，电感电流）为 ($1V$，$0A$)，但在虚平面上应表示成($1V$，$iA$)，其实部为($1V$，$0A$)，后续将进一步揭示这种表现形式的原因。
若选择合适的控制周期，使对应的离散系统状态转移矩阵$\varphi$的特征值为$i$，即$e^{\frac{\pi }{2}}$，那么经过一拍的时间后，系统状态变成了($iV$，$-1A$)，取其实数部分，得物理状态($0V$，$-1A$)，说明电容的能量全部转移到电感上了。再过一拍，系统状态为($-1V$,$-iA$), 取实部得物理状态($-1V$, $0A$)。这样随着时间推移，接下来的系统状态为 ($-iV$, $1A$），($1V$, $iA$) 等等，其实部为实际物理量，如此往复循环。

如果我们把控制周期缩短，则特征值将会改变。若特征值为$e^{\frac{\pi }{4}}$，那么我们就可以更清楚地看到系统状态变化的过程如下表所示。可以看出，电容电压、电感电流都如正弦振荡般往复循环。

无论在复频域中或者是Z域中，但凡特征值为复数，都意味着系统有一对共轭极点。共轭极点的位置表示系统的固有频率及其幅度变化情况。为方便理解，上面以Z平面中的特点举例（采样时间T不同，Z变换的结果也不同，那么特征值也不一样）。若特征值都在单位圆上，因此系统输出会等幅振荡，这是上述的情况。如果特征值在单位圆外，振荡幅值将不断增加（比如自激振荡之类），特征值在单位圆内，振荡幅值会不断衰减（把上述电路串一个电阻）。

上面讨论的是系统状态与特征向量方向相同的情况。但因为系统状态可以拆分成多个特征向量的线性组合，因而该论述对于任何的系统初始状态，都是适用的。

最后，给大家推荐一本《线性代数的几何意义》，是它带给了我灵感。