矩阵运算

#include <iostream>
#include <cassert>
#include <cmath>
#include <vector>
#include <string.h>

#define MAX 20
#define E 0.000000001

using namespace std;

// 定义向量别名
typedef vector<double> Vec;

//有关向量运算需要重载相关运算符，双目操作符左右值若是矩阵则调用，不是矩阵则调用非重载版本

// 重载-运算符
Vec operator-(const Vec& x, const Vec& y)
{
assert(x.size() == y.size());
// #include <cassert>
Vec tmp;
for (size_t i = 0; i < x.size(); ++i)
tmp.push_back(x[i] - y[i]);
return tmp;         // 返回局部变量的拷贝
}

// 重载*运算符  求内积
double operator*(const Vec& x, const Vec& y)
{
assert(x.size() == y.size());       // #include <cassert>
double sum = 0.;
for (size_t i = 0; i < x.size(); ++i)
sum += x[i] * y[i];
return sum;
}

// 求外积
// 三维的情况
Vec operator^(const Vec& x, const Vec& y)
{
assert(x.size() == y.size() && x.size() == 3);
return Vec{ x[1] * y[2] - x[2] * y[1],
x[2] * y[0] - x[0] * y[2],
x[0] * y[1] - x[1] * y[0] };
// uniform initialization, C++11新特性（这里放在VC6.0不可运行，注意一下）
}

// 二维就姑且返回其模长吧
double twoDCrossProd(const Vec& x, const Vec& y)
{
return x[0] * y[1] - x[1] * y[0];
}

// 模长或者范数的计算
double norm(const Vec& x)
{
double val = 0.;
for (auto elem : x)
val += elem*elem;
return sqrt(val);
// #include <cmath>
}

//求矩阵的转置
void swap(double *a, double *b)
{
double tmp = *a;
*b = *a;
*a = tmp;
}

void matrix_transpose(vector<Vec>& vec)
{
int i, j;
int n = vec[0].size();
for (i = 1; i<n; i++)
{
for (j = 0; j<i; j++)
swap(vec[i][j], vec[j][i]);
}
}

//从vector转到C形式的矩阵
void vector_c(double dst[MAX][MAX], const vector<Vec> &vec, int n)
{
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++)
{
dst[i][j] = vec[i][j];
}
}

/**
* 计算矩阵src的模
*/
double calculate_A(double src[][MAX], int n)
{
int i, j, k, x, y;
double tmp[MAX][MAX], t;
double result = 0.0;

if (n == 1)
{
return src[0][0];
}

for (i = 0; i < n; ++i)
{
for (j = 0; j < n - 1; ++j)
{
for (k = 0; k < n - 1; ++k)
{
x = j + 1;
y = k >= i ? k + 1 : k;

tmp[j][k] = src[x][y];
}
}

t = calculate_A(tmp, n - 1);

if (i % 2 == 0)
{
result += src[0][i] * t;
}
else
{
result -= src[0][i] * t;
}
}

return result;
}

/**
* 计算伴随矩阵
*/
void calculate_A_adjoint(double src[MAX][MAX], double dst[MAX][MAX], int n)
{
int i, j, k, t, x, y;
double tmp[MAX][MAX];

if (n == 1)
{
dst[0][0] = 1;
return;
}

for (i = 0; i < n; ++i)
{
for (j = 0; j < n; ++j)
{
for (k = 0; k < n - 1; ++k)
{
for (t = 0; t < n - 1; ++t)
{
x = k >= i ? k + 1 : k;
y = t >= j ? t + 1 : t;

tmp[k][t] = src[x][y];
}
}

dst[j][i] = calculate_A(tmp, n - 1);

if ((i + j) % 2 == 1)
{
dst[j][i] = -1 * dst[j][i];
}
}
}
}

/**
* 得到逆矩阵
*/
bool calculate_A_inverse(double src[MAX][MAX], double dst[MAX][MAX], int n)
{
double A = calculate_A(src, n);
double tmp[MAX][MAX];
int i, j;

if (fabs(A - 0) <= E)
{
printf("不可能有逆矩阵！\n");
return false;
}

calculate_A_adjoint(src, tmp, n);

for (i = 0; i < n; ++i)
{
for (j = 0; j < n; ++j)
{
dst[i][j] = (double)(tmp[i][j] / A);
}
}

return true;
}

/**
* 输出矩形查看
*/
void print_M(double M[][MAX], int n)
{
int i, j;

for (i = 0; i < n; ++i)
{
for (j = 0; j < n; ++j)
{
printf("%.2f ", M[i][j]);
}

printf("\n");
}
}

//输出一维向量
void print_vector(const Vec &vec)
{
int len = vec.size();
for (int i = 0; i < len - 1; i++)
cout << vec[i] << ' ';
cout << vec[len - 1] << endl;
}

//输出二维向量
void print_vector2(const vector<Vec>& vec)
{
int i, j;
int n = vec[0].size();
for (i = 0; i < n; ++i)
{
for (j = 0; j < n; ++j)
{
printf("%.2f ", vec[i][j]);
}

printf("\n");
}
}

//求*积，两矩阵相乘
void calculate_product(double product[MAX][MAX], double A[MAX][MAX], double B[MAX][MAX], int n)
{
for (int i = 0; i<n; i++)
for (int j = 0; j<n; j++)
for (int k = 0; k<n; k++)
product[i][j] += A[i][k] * B[k][j];

}

/**
* main
*/
int main()
{
/**
* 测试矩阵
*/

cout << "****** 终极目标:R=M*S ********" << endl;
cout << endl;

////////////////////////////////////////////////////////////////////////
//创建M

cout << "####首先求M矩阵####" << endl;

cout << "*****M矩阵相关的三个一维向量****" << endl;
//输入两个一维向量P1=(x1,y1,z1),P2=(x2,y2,z2),P3=(x3,y3,z3) N为列数
int M_N;
cout << "请输入三个一维向量列数:";
cin >> M_N;

//创建三个一维向量
Vec P1(M_N, 0);
Vec P2(M_N, 0);
Vec P3(M_N, 0);

cout << "请输入第一个一维向量:";
for (int i = 0; i < M_N; i++)
cin >> P1[i];
cout << "请输入第二个一维向量:";
for (int i = 0; i < M_N; i++)
cin >> P2[i];
cout << "请输入第三个一维向量:";
for (int i = 0; i < M_N; i++)
cin >> P3[i];

cout << "\n****M成功创建****" << endl;

//求两个向量差
cout << "P2P1 = P1 - P2：" << endl;
Vec P2P1 = P1 - P2;
print_vector(P2P1);

cout << "P2P3 = P3 - P2：" << endl;
Vec P2P3 = P3 - P2;
print_vector(P2P3);

//求P2P1和P2P3两个向量的外积
cout << "Outer_product = P2P1 X P2P3：" << endl;
Vec Outer_product_M = P2P1^P2P3;
print_vector(Outer_product_M);

//构造M矩阵
cout << "****构造M矩阵****" << endl;
vector<Vec> src_vector_M{ P2P1, P2P3, Outer_product_M };//C++11新特性
//vector形式
cout << "vector形式：" << endl;
matrix_transpose(src_vector_M);
print_vector2(src_vector_M);

//dst[][]形式
double dst_M[MAX][MAX];
cout << "dst[][]形式：" << endl;
int len = src_vector_M[0].size();
vector_c(dst_M, src_vector_M, len);
print_M(dst_M, len);

cout << "****构造M矩阵结束****" << endl;
cout << endl;

////////////////////////////////////////////////////////////////////////
//创建N

cout << "####然后求S矩阵####" << endl;
cout << "*****N矩阵相关的三个一维向量****" << endl;
//输入两个一维向量Q1=(x1,y1,z1),Q2=(x2,y2,z2),Q3=(x3,y3,z3) N为列数
int N_N;
cout << "请输入三个一维向量列数:";
cin >> N_N;

//创建三个一维向量
Vec Q1(N_N, 0);
Vec Q2(N_N, 0);
Vec Q3(N_N, 0);

cout << "请输入第一个一维向量:";
for (int i = 0; i < N_N; i++)
cin >> Q1[i];
cout << "请输入第二个一维向量:";
for (int i = 0; i < N_N; i++)
cin >> Q2[i];
cout << "请输入第三个一维向量:";
for (int i = 0; i < N_N; i++)
cin >> Q3[i];

cout << "\n****N成功创建****" << endl;

//求两个向量差
cout << "Q2Q1 = Q1 - Q2：" << endl;
Vec Q2Q1 = Q1 - Q2;
print_vector(Q2Q1);

cout << "Q2Q3 = Q3 - Q2：" << endl;
Vec Q2Q3 = Q3 - Q2;
print_vector(Q2Q3);

//求P2P1和P2P3两个向量的外积
cout << "Outer_product_N = Q2Q1 X Q2Q3：" << endl;
Vec Outer_product_N = Q2Q1^Q2Q3;
print_vector(Outer_product_N);

//构造N矩阵
cout << "****构造N矩阵****" << endl;
vector<Vec> src_vector_N{ Q2Q1, Q2Q3, Outer_product_N };//C++11新特性
//vector形式
cout << "vector形式：" << endl;
matrix_transpose(src_vector_N);
print_vector2(src_vector_N);

//dst[][]形式
double dst_N[MAX][MAX];
cout << "dst[][]形式：" << endl;
int length = src_vector_N[0].size();
vector_c(dst_N, src_vector_N, length);
print_M(dst_N, length);

cout << "****构造N矩阵结束****" << endl;

//求M的逆矩阵S
double S[MAX][MAX];
bool is_exist = calculate_A_inverse(S, dst_M, length);

if (is_exist)
{
print_M(S, length);
}
else
{
printf("不存在!\n");
}

////////////////////////////////////////////////////////////////////////
////求终极目标R
cout << "####求终极目标R####" << endl;
double R[MAX][MAX];
calculate_product(R, dst_M, S, length);
print_M(R, length);

///////////逆矩阵函数测试区，把无关代码屏蔽/////////////////
//double test[MAX][MAX], dst[MAX][MAX];
//int n = 3;

///**
//* 构造最简单的：
//*  1, 0, 0,
//*  0, 2, 0,
//*  0, 0, 5
//*/
//memset(test, 0, sizeof(test));

//test[0][0] = 1.0;
//test[1][1] = 2.0;
//test[2][2] = 5.0;

//bool is_exist = calculate_A_inverse(test, dst, n);

//if (is_exist)
//{
//	print_M(dst, n);
//}
//else
//{
//	printf("不存在!\n");
//}

system("pause");
return 0;
}