行列式求值、矩阵求逆

#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 rowsfor (col = 0; col < n; col++){result = mat[i][col];mat[i][col] = mat[j][col];mat[j][col] = result;}sign = -sign;}//sub i rowfor (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 << "输入行列式的阶数：";}
}