矩阵运算

#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;
}