#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 << "输入行列式的阶数:";
}
}
行列式求值、矩阵求逆
最新推荐文章于 2024-02-04 21:42:24 发布