#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<<"输入行列式的阶数:";
}
}
矩阵操作
最新推荐文章于 2022-01-18 17:38:33 发布