行列式求值、矩阵求逆_cuda矩阵求逆

#include <iostream>
#include <string>
#include <assert.h>
#include <malloc.h>
#include <iostream>
#include <stdlib.h>
#include <memory.h>
#include <time.h>
 
using namespace std;
 
//动态分配大小位size的一维数组
template<typename T>
bool allocateMemory1D(T **p, const int size)
{
	*p = NULL;
	*p = (T *)malloc(size * sizeof(T));
 
	if (*p == NULL) return false;
	return true;
}
 
//动态分配rows行、cols列的二维数组
template<typename T>
bool allocateMemory2D(T ***p, const int rows, const int cols)
{
	*p = NULL;
	*p = (T **)malloc(rows * sizeof(T *));
	if (*p == NULL) return false;
 
	memset(*p, 0, sizeof(T *)* rows);
	for (int i = 0; i < rows; i++)
	{
		(*p)[i] = (T *)malloc(cols * sizeof(T));
		if ((*p)[i] == NULL) return false;
	}
 
	return true;
}
 
//释放一维动态数组的空间
template<typename T>
void freeMemory1D(T **p)
{
	if (*p != NULL)
	{
		free(*p);
		*p = NULL;
	}
}
 
//释放二维动态数组的空间
template<typename T>
void freeMemory2D(T ***p, const int rows)
{
	if (*p != NULL)
	{
		for (int i = 0; i < rows; i++)
		{
			if ((*p)[i] != NULL)
			{
				free((*p)[i]);
				(*p)[i] = NULL;
			}
		}
 
		free(*p);
		*p = NULL;
	}
}
 
//template<class T>
//void swap(T &a, T &b)
//{
//	T tmp = a;
//	a = b;
//	b = tmp;
//}
 
//用随机数填充二维数组(矩阵)
template<typename T, typename V>
bool generateNNumbers2D(T **a, int rows, int cols, V minValue, V maxValue)
{
	//assert(a != NULL && rows > 0 && cols > 0);
	if (a == NULL || rows < 1 && cols < 1) return false;
	if (minValue > maxValue) swap(minValue, maxValue);
	const int N = 999999;
 
	srand((unsigned)time(NULL));
	for (int i = 0; i < rows; i++)
	{
		for (int j = 0; j < cols; j++)
		{
			a[i][j] = rand() % (int)(maxValue - minValue + 1) + minValue + rand() % (N + 1) / (float)(N + 1);
		}
	}
 
	return true;
}
 
void displayArray2D(double **a, const int rows, const int cols)
{
	assert(a != NULL && rows > 0 && cols > 0);
 
	for (int i = 0; i < rows; i++)
	{
		printf("%d: ", i + 1);
		for (int j = 0; j < cols; j++)
		{
			printf("%lf ", a[i][j]);
		}
		printf("\n");
	}
}
 
//降阶法递归求行列式的值，就是按照线性代数书上的公式，我是按照第一行进行展开
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 = NULL;
 
		allocateMemory2D(&tmpMat, n - 1, n - 1);
 
		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;
		}
 
		freeMemory2D(&tmpMat, n - 1);
 
		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;
}
 
//求一个矩阵的邻接矩阵,T为输入矩阵，adjointMat存放T对应的邻接矩阵，n为矩阵的阶数
template <typename T>
bool adjointMatrix(T **mat, T ** &adjointMat, int n)
{
	int i, j;
	if (mat == NULL || n < 1) return false;
 
	if (n == 1)
	{
		adjointMat[0][0] = 1;
		return true;
	}
 
	T **tmpMat = NULL;
	allocateMemory2D(&tmpMat, n - 1, 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;
		}
	}
 
	freeMemory2D(&tmpMat, n - 1);
 
	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;
}
 
//两个矩阵相乘，mat3 = mat1 * mat2
template <typename T>
bool matrixMultiply(T **mat1, T **mat2, T **&mat3, int n)
{
	if (mat1 == NULL || mat2 == NULL || n < 1) return false;
	
	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];
		}
	}
 
	return true;
}
 
void main(void)
{
	int n, i, j, ret;
	double **mat = NULL;
	double **adjointMat = NULL;
	double absouluteError, relativeError, d;
	double **tmpMat = NULL;
	bool isReversible;
	bool isDataAutoGenerate;
	char ch;
	char buf[256];
 
	cout << "自动生成矩阵吗(y/Y：是)：";
	cin >> ch;
	isDataAutoGenerate = (ch == 'y' || ch == 'Y');
	cin.getline(buf, sizeof(buf));	//清空输入缓冲区
 
	cout << "输入行列式的阶数：";
	while (cin >> n)
	{
		allocateMemory2D(&mat, n, n);
		if (isDataAutoGenerate)
		{
			generateNNumbers2D(mat, n, n, -100, 100);
			cout << "自动生成的" << n << "阶矩阵为：" << endl;
			displayArray2D(mat, n, n);
 
		}
		else
		{
			cout << "依次输入矩阵的每一个元素：\n";
			for (i = 0; i < n * n; i++) cin >> mat[i / n][i % n];
		}
 
		cout << "分别用两种方法计算行列式的值."<<endl;
		d = det(mat, n);
		cout << "1.行列式的值为：" << d << endl;
		cout << "2.行列式的值为：" << det1(mat, n) << endl;
 
		allocateMemory2D(&adjointMat, n, n);
		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;
			}
		}
 
		ret = inverseMatrix(d, adjointMat, adjointMat, n);
		if (ret == -2) 
		{
			isReversible = false;
			cout << "参数错误" << endl;
		}
		else if (ret == -1)
		{
			isReversible = false;
			cout << "矩阵不可逆" << endl;
		}
		else
		{
			isReversible = true;
			cout << "逆矩阵:" << endl;
 
			for (i = 0; i < n; i++)
			{
				for (j = 0; j < n; j++)
				{
					cout<<adjointMat[i][j]<<" ";
				}
				cout<<endl;
			}
		}
 
		if (isReversible)
		{
			cout << "计算一个矩阵与其对应的逆矩阵的乘积：" << endl;
			allocateMemory2D(&tmpMat, n, n);
			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 (isReversible)
		{
			cout << "原矩阵 * 逆矩阵 和 原矩阵同规模的单位矩阵作对比，看误差多大." << endl;
			compareMatrix(mat, tmpMat, n, &absouluteError, &relativeError);
			cout << "绝对误差：" << absouluteError << endl;
			cout << "相对误差：" << relativeError << endl;
		}
 
		freeMemory2D(&mat, n);
		freeMemory2D(&adjointMat, n);
		freeMemory2D(&tmpMat, n);
 
		cout << "输入行列式的阶数：";
	}
}