矩阵操作

#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 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;}//求一个矩阵的邻接矩阵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<<"输入行列式的阶数：";}}