C語言通過QR分解計算矩陣的特征值和特征向量

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

//--------------------------這裡是一些定義的結構體和資料類型---------
//純C裡面定義的布爾類型
typedef enum { False = 0, True = 1 }Bool;
//定義矩陣元素的類型為matrixType
typedef double matrixType;

//此結構體用來表示矩陣，其中row為行，column為列，height為高，array用來存放矩陣元素(用一維來模拟二維/三維)
typedef struct
{
	unsigned  row, column, height;
	matrixType *array; //使用時，必須對*array進行初始化
}Matrix;

//---------下面是QR分解，求解線性方程所用到的一些函數-----------
/*
matrix為要設定大小并配置設定記憶體的矩陣，row、column、height分别為行，列，高。
函數調用成功則則傳回true,否則傳回false
*/
Bool SetMatrixSize(Matrix *matrix, const unsigned row, const unsigned column, const unsigned height)
{
	unsigned size = row * column * height * sizeof(matrixType);
	if (size <= 0)
	{
		return False;
	}

	matrix->array = (matrixType*)malloc(size);

	//如果配置設定記憶體成功
	if (matrix->array)
	{
		matrix->row = row;
		matrix->column = column;
		matrix->height = height;
		return True;
	}
	else
	{
		matrix->row = matrix->column = matrix->height = 0;
		return False;
	}
}

//設定Matrix矩陣中的所有元素大小為ele
void SetMatrixEle(Matrix *matrix, matrixType ele)
{
	unsigned size = matrix->row * matrix->column * matrix->height;
	unsigned i;

	for (i = 0; i < size; ++i)
	{
		matrix->array[i] = ele;
	}
}

//設定Matrix矩陣中的所有元素大小為0
void SetMatrixZero(Matrix*matrix)
{
	SetMatrixEle(matrix, 0);
}

//判斷矩陣是否為空，若為空則傳回1，否則傳回0
Bool IsNullMatrix(const Matrix* matrix)
{
	unsigned size = matrix->row * matrix->column * matrix->column;

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

//銷毀矩陣，即釋放為矩陣動态配置設定的記憶體,并将矩陣的長寬高置0
void DestroyMatrix(Matrix *matrix)
{
	if (!IsNullMatrix(matrix))
	{
		free(matrix->array);
		matrix->array = NULL;
	}

	matrix->row = matrix->column = matrix->height = 0;
}

//計算矩陣可容納元素個數，即return row*column*height
unsigned MatrixCapacity(const Matrix*matrix)
{
	return matrix->row * matrix->column * matrix->height;
}


//||matrix||_2  求A矩陣的2範數
matrixType MatrixNorm2(const Matrix *matrix)
{
	unsigned size = matrix->row * matrix->column *matrix->height;
	unsigned i;
	matrixType norm = 0;

	for (i = 0; i < size; ++i)
	{
		norm += (matrix->array[i]) *(matrix->array[i]);
	}

	return (matrixType)sqrt(norm);
}

//matrixB = matrix(:,:,height)即拷貝三維矩陣的某層，若matrix為二維矩陣，需将height設定為0
Bool CopyMatrix(Matrix* matrixB, Matrix *matrix, unsigned height)
{
	unsigned size, i;
	Matrix matrixA;

	//判斷height值是否正确
	if (height < 0 || height >= matrix->height)
	{
		printf("ERROR: CopyMatrix() parameter error！\n");
		return False;
	}

	//将matrix(:,:,height) 轉換為二維矩陣matrixA
	matrixA.row = matrix->row;
	matrixA.column = matrix->column;
	matrixA.height = 1;
	matrixA.array = matrix->array + height * matrix->row * matrix->column;

	//判斷兩矩陣指向的記憶體是否相等
	if (matrixA.array == matrixB->array)
	{
		return True;
	}

	//計算matrixA的容量
	size = MatrixCapacity(&matrixA);
	//判斷matrixB的容量與matrixA的容量是否相等
	if (MatrixCapacity(matrixB) != size)
	{
		DestroyMatrix(matrixB);
		SetMatrixSize(matrixB, matrixA.row, matrixA.column, matrixA.height);
	}
	else
	{
		matrixB->row = matrixA.row;
		matrixB->column = matrixA.column;
		matrixB->height = matrixA.height;
	}

	for (i = 0; i < size; ++i)
	{
		matrixB->array[i] = matrixA.array[i];
	}

	return True;
}

//matrixC = matrixA * matrixB
Bool MatrixMulMatrix(Matrix *matrixC, const Matrix* matrixA, const Matrix* matrixB)
{
	size_t row_i, column_i, i;
	size_t indexA, indexB, indexC;
	matrixType temp;
	Matrix tempC;

	if (IsNullMatrix(matrixA) || IsNullMatrix(matrixB))
	{
		return False;
	}

	if (matrixA->column != matrixB->row)
	{
		return False;
	}

	if (MatrixCapacity(matrixC) != matrixA->row * matrixB->column)
	{
		SetMatrixSize(&tempC, matrixA->row, matrixB->column, 1);
	}
	else
	{
		tempC.array = matrixC->array;
		tempC.row = matrixA->row;
		tempC.column = matrixB->column;
		tempC.height = 1;
	}

	for (row_i = 0; row_i < tempC.row; ++row_i)
	{
		for (column_i = 0; column_i < tempC.column; ++column_i)
		{
			temp = 0;
			for (i = 0; i < matrixA->column; ++i)
			{
				indexA = row_i * matrixA->column + i;
				indexB = i * matrixB->column + column_i;

				temp += matrixA->array[indexA] * matrixB->array[indexB];
			}
			indexC = row_i * tempC.column + column_i;

			tempC.array[indexC] = temp;
		}
	}


	if (tempC.array != matrixC->array)
	{
		DestroyMatrix(matrixC);

		matrixC->array = tempC.array;
	}

	matrixC->row = tempC.row;
	matrixC->column = tempC.column;
	matrixC->height = tempC.height;



	return True;
}

//對vector中所有元素排序，sign= 0 時為升序，其餘為降序
void SortVector(Matrix *vector, int sign)
{
	matrixType mid;
	int midIndex;
	int size = MatrixCapacity(vector);
	int i, j;

	if (0 == sign)
	{
		for (i = 0; i < size; ++i)
		{
			mid = vector->array[i];
			midIndex = i;
			for (j = i + 1; j < size; ++j)
			{
				if (mid > vector->array[j])
				{
					mid = vector->array[j];
					midIndex = j;
				}
			}

			vector->array[midIndex] = vector->array[i];
			vector->array[i] = mid;
		}
	}
	else
	{
		for (i = 0; i < size; ++i)
		{
			mid = vector->array[i];
			midIndex = i;
			for (j = i + 1; j < size; ++j)
			{
				if (mid < vector->array[j])
				{
					mid = vector->array[j];
					midIndex = j;
				}
			}

			vector->array[midIndex] = vector->array[i];
			vector->array[i] = mid;
		}
	}
}

//列印矩陣
void PrintMatrix(const Matrix *matrix)
{
	size_t row_i, column_i, height_i, index;

	for (height_i = 0; height_i < matrix->height; ++height_i)
	{
		(matrix->height == 1) ? printf("[:,:] = \n") : printf("[%d,:,:] = \n", height_i);

		for (row_i = 0; row_i < matrix->row; ++row_i)
		{
			for (column_i = 0; column_i < matrix->column; ++column_i)
			{
				index = height_i * matrix->row * matrix->column + row_i * matrix->column + column_i;
				printf("%12.4g", matrix->array[index]);
			}
			printf("\n");
		}
	}
}

//----------------------QR分解-------------------------------------------

//将A分解為Q和R
void QR(Matrix *A, Matrix *Q, Matrix *R)
{
	unsigned  i, j, k, m;
	unsigned size;
	const unsigned N = A->row;
	matrixType temp;

	Matrix a, b;

	//如果A不是一個二維方陣，則提示錯誤，函數計算結束
	if (A->row != A->column || 1 != A->height)
	{
		printf("ERROE: QR() parameter A is not a square matrix!\n");
		return;
	}

	size = MatrixCapacity(A);
	if (MatrixCapacity(Q) != size)
	{
		DestroyMatrix(Q);
		SetMatrixSize(Q, A->row, A->column, A->height);
		SetMatrixZero(Q);
	}
	else
	{
		Q->row = A->row;
		Q->column = A->column;
		Q->height = A->height;
	}

	if (MatrixCapacity(R) != size)
	{
		DestroyMatrix(R);
		SetMatrixSize(R, A->row, A->column, A->height);
		SetMatrixZero(R);
	}
	else
	{
		R->row = A->row;
		R->column = A->column;
		R->height = A->height;
	}

	SetMatrixSize(&a, N, 1, 1);
	SetMatrixSize(&b, N, 1, 1);

	for (j = 0; j < N; ++j)
	{
		for (i = 0; i < N; ++i)
		{
			a.array[i] = b.array[i] = A->array[i * A->column + j];
		}

		for (k = 0; k < j; ++k)
		{
			R->array[k * R->column + j] = 0;

			for (m = 0; m < N; ++m)
			{
				R->array[k * R->column + j] += a.array[m] * Q->array[m * Q->column + k];
			}

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

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

		for (i = 0; i < N; ++i)
		{
			Q->array[i * Q->column + j] = b.array[i] / temp;
		}
	}

	DestroyMatrix(&a);
	DestroyMatrix(&b);
}

//----------------------使用特征值計算矩陣特征向量-----------------
//eigenVector為計算結果即矩陣A的特征向量
//eigenValue為矩陣A的所有特征值，
//A為要計算特征向量的矩陣
void Eigenvectors(Matrix *eigenVector, Matrix *A, Matrix *eigenValue)
{
	unsigned i, j, q;
	unsigned count;
	int m;
	const unsigned NUM = A->column;
	matrixType eValue;
	matrixType sum, midSum, mid;
	Matrix temp;

	SetMatrixSize(&temp, A->row, A->column, A->height);

	for (count = 0; count < NUM; ++count)
	{
		//計算特征值為eValue，求解特征向量時的系數矩陣
		eValue = eigenValue->array[count];
		CopyMatrix(&temp, A, 0);
		for (i = 0; i < temp.column; ++i)
		{
			temp.array[i * temp.column + i] -= eValue;
		}

		//将temp化為階梯型矩陣
		for (i = 0; i < temp.row - 1; ++i)
		{
			mid = temp.array[i * temp.column + i];
			for (j = i; j < temp.column; ++j)
			{
				temp.array[i * temp.column + j] /= mid;
			}

			for (j = i + 1; j < temp.row; ++j)
			{
				mid = temp.array[j * temp.column + i];
				for (q = i; q < temp.column; ++q)
				{
					temp.array[j * temp.column + q] -= mid * temp.array[i * temp.column + q];
				}
			}
		}
		midSum = eigenVector->array[(eigenVector->row - 1) * eigenVector->column + count] = 1;
		for (m = temp.row - 2; m >= 0; --m)
		{
			sum = 0;
			for (j = m + 1; j < temp.column; ++j)
			{
				sum += temp.array[m * temp.column + j] * eigenVector->array[j * eigenVector->column + count];
			}
			sum = -sum / temp.array[m * temp.column + m];
			midSum += sum * sum;
			eigenVector->array[m * eigenVector->column + count] = sum;
		}

		midSum = sqrt(midSum);
		for (i = 0; i < eigenVector->row; ++i)
		{
			eigenVector->array[i * eigenVector->column + count] /= midSum;
		}
	}
	DestroyMatrix(&temp);
}
int main()
{
	const unsigned NUM = 50; //最大疊代次數

	unsigned N = 4;
	unsigned k;

	Matrix A, Q, R, temp;
	Matrix eValue;


	//配置設定記憶體
	SetMatrixSize(&A, N, N, 1);
	SetMatrixSize(&Q, A.row, A.column, A.height);
	SetMatrixSize(&R, A.row, A.column, A.height);
	SetMatrixSize(&temp, A.row, A.column, A.height);
	SetMatrixSize(&eValue, A.row, 1, 1);

	//A設定為一個簡單矩陣
	A.array[0] = -1;
	A.array[1] = 2;
	A.array[2] = 1;
	A.array[3] = 2;
	A.array[4] = 4;
	A.array[5] = -1;
	A.array[6] = 1;
	A.array[7] = 1;
	A.array[8] = -6;
	A.array[9] = -6;
	A.array[10] = -6;
	A.array[11] = -6;
	A.array[12] = -6;
	A.array[13] = -6;
	A.array[14] = -6;
	A.array[15] = -6;


	//拷貝A矩陣元素至temp
	CopyMatrix(&temp, &A, 0);

	//初始化Q、R所有元素為0
	SetMatrixZero(&Q);
	SetMatrixZero(&R);
	//使用QR分解求矩陣特征值
	for (k = 0; k < NUM; ++k)
	{
		QR(&temp, &Q, &R);
		MatrixMulMatrix(&temp, &R, &Q);
	}
	//擷取特征值，将之存儲于eValue
	for (k = 0; k < temp.column; ++k)
	{
		eValue.array[k] = temp.array[k * temp.column + k];
	}

	//對特征值按照降序排序
	SortVector(&eValue, 1);

	//根據特征值eValue，原始矩陣求解矩陣特征向量Q
	Eigenvectors(&Q, &A, &eValue);

	//列印特征值
	printf("特征值：");
	PrintMatrix(&eValue);

	//列印特征向量
	printf("特征向量");
	PrintMatrix(&Q);
	DestroyMatrix(&A);
	DestroyMatrix(&R);
	DestroyMatrix(&Q);
	DestroyMatrix(&eValue);
	DestroyMatrix(&temp);

	return 0;
}