矩阵操作

#include <iostream>
#include <string>
#include <assert.h>
#include <malloc.h>
#include <iostream>
#include <stdlib.h>
#include <memory.h>

using namespace std;

//降阶法求行列式的值，就是按照线性代数书上的公式，我是按照第一行进行展开
template <typename T>
double static det(T **mat, const int n)
{
	assert(mat != NULL && n > 0);

	if (n == 1) return (double)mat[0][0];
	else if (n == 2) return (double)(mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0]);
	else 
	{
		int i, j, k, flag = 1, col;
		double value = 0.0;

		T **tmpMat = (T **)malloc(sizeof(T *) * (n - 1));
		assert(tmpMat != NULL);
		memset(tmpMat, 0, sizeof(T *) * (n - 1));

		for (i = 0; i < n - 1; i++)
		{
			tmpMat[i] = (T *)malloc(sizeof(T) * (n - 1));
			assert(tmpMat[i] != NULL);
		}

		for (i = 0; i < n; i++)
		{
			for (j = 1; j < n; j++)
			{
				col = 0;
				for (k = 0; k < n; k++)
				{
					if (k != i)
					{
						tmpMat[j - 1][col++] = mat[j][k];
					}
				}
			}

			value += mat[0][i] * det(tmpMat, n - 1) * flag;
			flag = -flag;
		}

		for (i = 0; i < n - 1; i++)
		{
			if(tmpMat[i] != NULL) free(tmpMat[i]);
		}
		if(tmpMat != NULL) free(tmpMat);
		
		return value;
	}
}

//将转换为上三角形式后再求行列式的值
double static det1(double **mat, const int n)
{
	assert(mat && n > 0);
	const double PRECESION = 1E-6;
	int row, col, i, j;
	bool flag = false;
	int sign = 1;
	double result = 1.0;

	for (i = 0; i < n - 1; i++)
	{
		for (j = i; j < n; j++)
		{
			if (fabs(mat[i][j]) > PRECESION) break;
		}

		if(j >= n)
		{
			flag = true;
			break;
		}

		if (j != i)
		{
			//swap rows
			for (col = 0; col < n; col++)
			{
				result = mat[i][col];
				mat[i][col] = mat[j][col];
				mat[j][col] = result;
			}
			sign = -sign;
		}

			//sub i row
		for (row = j + 1; row < n; row++)
		{
			double base = mat[row][i] / mat[i][i];
			for (col = 0; col < n; col++)
			{
				mat[row][col] -= mat[i][col] * base;
			}
		}
	}

	if (flag)
	{
		return 0;
	}

	else 
	{
		result = 1.0;
		for (i = 0; i < n; i++)
		{
			result *= mat[i][i];
		}
		if (sign < 0)
		{
			result = -result;
		}
	}
	
	return result;
}

//求一个矩阵的邻接矩阵
template <typename T>
bool adjointMatrix(T **mat, T ** &adjointMat, int n)
{
	int i, j;
	if (mat == NULL || n < 1) return false;

	if (adjointMat == NULL)
	{
		adjointMat = new T*[n];
		for (i = 0; i < n; i++)
		{
			adjointMat[i] = new T[n];
		}
	}

	if (n == 1)
	{
		adjointMat[0][0] = 1;
		return true;
	}

	T **tmpMat = NULL;
	tmpMat = new T*[n - 1];
	memset(tmpMat, 0, sizeof(*tmpMat) * (n - 1));
	for (i = 0; i < n - 1; i++)
	{
		tmpMat[i] = new T[n - 1];
	}

	int sign = -1, row, col, rowIndex, colIndex;
	for (i = 0; i < n; i++)
	{
		sign = -sign;
		int s = sign;
		for (j = 0; j < n; j++)
		{
			rowIndex = 0;

			for (row = 0; row < n; row++)
			{
				colIndex = 0;

				if (row != i)
				{
					for (col = 0; col < n; col++)
					{
						if (col != j)
						{
							tmpMat[rowIndex][colIndex] = mat[row][col];
							
							colIndex++;
						}
					}

					rowIndex++;
				}
			}

			adjointMat[j][i] = s * det(tmpMat, n - 1);
			
			s = -s;
		}
	}

	for (i = 0; i < n - 1; i++) if(tmpMat[i]) delete[] tmpMat[i];
	if (tmpMat) delete[] tmpMat;

	return true;
}

//求一个矩阵的逆矩阵
template <typename T>
int inverseMatrix(double d, T **mat, T ** &inverseMat, int n)
{
	if (n < 1 || mat == NULL || inverseMat == NULL) return -2;

	//double d = det(mat, n);
	if (fabs(d) < 1E-6) return -1;

	for (int row = 0; row < n; row++)
	{
		for (int col = 0; col < n; col++)
		{
			inverseMat[row][col] = mat[row][col] / d;
		}
	}

	return 0;
}

//两个矩阵相乘
template <typename T>
bool matrixMultiply(T **mat1, T **mat2, T **&mat3, int n)
{
	if (mat1 == NULL || mat2 == NULL || n < 1) return false;

	if(mat3 == NULL)
	{
		mat3 = new T*[n];
		memset(mat3, 0, sizeof(*mat3) * n);

		for (int i = 0; i < n; i++) mat3[i] = new T[n];
	}
	
	for (int row = 0; row < n; row++)
	{
		for (int col = 0; col < n; col++)
		{
			double sum = 0.0;
			for (int k = 0; k < n; k++)
			{
				sum += mat1[row][k] * mat2[k][col];
			}

			mat3[row][col] = sum;
		}
	}

	return true;
}

template <typename T>
bool compareMatrix(T **referenceMat, T **mat, int n, double *absoluteError, double *relativeError)
{
	if(referenceMat == NULL || mat == NULL || n < 1 || absoluteError == NULL || relativeError == NULL) return false;

	*absoluteError = 0;
	*relativeError = 0;

	for (int row = 0; row < n; row++)
	{
		for (int col = 0; col < n; col++)
		{
			double tmp = fabs(mat[row][col] - referenceMat[row][col]);
			*absoluteError += tmp;

			if (fabs(referenceMat[row][col]) > 1E-6) *relativeError += tmp / referenceMat[row][col];
		}
	}
}

void main()
{
	int n, i, j, ret;
	double **mat = NULL;
	double **adjointMat = NULL;
	double absouluteError, relativeError, d;
	double **tmpMat = NULL;

	cout<<"输入行列式的阶数：";
	while (cin>>n)
	{
		mat = (double **)malloc(sizeof(double *) * n);
		assert(mat);
		memset(mat, 0, sizeof(double *) * n);

		for (i = 0; i < n; i++)
		{
			mat[i] = (double *)malloc(sizeof(double) * n);
			assert(mat[i]);
		}

		for (i = 0; i < n * n; i++) cin>>mat[i / n][i % n];

		d = det(mat, n);
		cout<<"行列式的值为：";
		cout<<d<<", "<<det1(mat, n)<<endl;

		cout<<"伴随矩阵为："<<endl;
		if (!adjointMatrix(mat, adjointMat, n)) cout<<"伴随矩阵求解失败！"<<endl;
		else
		{
			for (i = 0; i < n; i++)
			{
				for (j = 0; j < n; j++)
				{
					cout<<adjointMat[i][j]<<" ";
				}
				cout<<endl;
			}
		}

		cout<<"逆矩阵:"<<endl;
		ret = inverseMatrix(d, adjointMat, adjointMat, n);
		if (ret == -2) cout<<"参数错误"<<endl;
		else if (ret == -1) cout<<"矩阵不可逆"<<endl;
		else{
			for (i = 0; i < n; i++)
			{
				for (j = 0; j < n; j++)
				{
					cout<<adjointMat[i][j]<<" ";
				}
				cout<<endl;
			}
		}

		cout<<"计算一个矩阵与其对应的逆矩阵的乘积："<<endl;
		ret = matrixMultiply(mat, adjointMat, tmpMat, n);
		for (i = 0; i < n; i++)
		{
			for (j = 0; j < n; j++)
			{
				cout<<tmpMat[i][j]<<" ";
			}
			cout<<endl;
		}

		//将mat矩阵置为单位矩阵
		for (i = 0; i < n; i++)
		{
			for (j = 0; j < n; j++)
			{
				if (i == j)
				{
					mat[i][j] = 1;
				}
				else mat[i][j] = 0;
			}
		}

		if (!ret) cout << "参数错误" << endl;
		else
		{
			compareMatrix(mat, tmpMat, n, &absouluteError, &relativeError);
			cout << "绝对误差：" << absouluteError << endl;
			cout << "相对误差：" << relativeError << endl;
		}

		for(i = 0; i < n; i++) if(mat[i] != NULL) free(mat[i]);
		if(mat != NULL) {free(mat); mat = NULL;}

		for (i = 0; i < n; i++) if (adjointMat[i] != NULL) free(adjointMat[i]);
		if (adjointMat != NULL) { free(adjointMat); adjointMat = NULL; }

		for (i = 0; i < n; i++) if (tmpMat[i] != NULL) free(tmpMat[i]);
		if (tmpMat != NULL) { free(tmpMat); tmpMat = NULL; }

		cout<<"输入行列式的阶数：";
	}
}