C++矩阵算法

// 矩阵计算------------------------------------------------------------------------------
	vector> MatrixInit0(int rows, int columns);	// 用0初始化矩阵
	vector> MatrixT(vector> vectMatrix, int rows, int columns);	//矩阵转置
	double MatrixDet(const vector> &vectMatrix, int rows);	// 矩阵行列式
	vector> MatrixInv(const vector> &vectMatrix, int rows);//矩阵求逆
	vector> MatrixMul(const vector> &vectMatrixA, const vector> &vectMatrixB,	//矩阵相乘
		int rows1, int columns1, int rows2, int columns2);
// 矩阵计算------------------------------------------------------------------------------	


// 用0初始化矩阵
vector> MatrixInit0(int rows, int columns)
{
	vector> vectMatrix;
	vectMatrix.resize(rows);
	for (int i = 0; i < rows; i++)
		vectMatrix[i].resize(columns, 0);
	return vectMatrix;
}

//矩阵转置
vector> MatrixT(vector> vectMatrix, int rows, int columns)
{
	vector> vectMatrixT;
	vectMatrixT = MatrixInit0(columns, rows);
	for (int i = 0; i < columns; i++)
		for (int j = 0; j < rows; j++)
			vectMatrixT[i][j] = vectMatrix[j][i];
	return vectMatrixT;
}

// 矩阵行列式
double MatrixDet(const vector> &vectMatrix, int rows)
{
	if (rows == 1)
		return vectMatrix[0][0];
	else {
		double dSum = 0;
		for (int j = 0; j < rows; j++) {
			vector> vectMatrixS;
			for (int i = 1; i < rows; i++) {
				vector vectTemp;
				for (int z = 0; z < rows; z++) {
					if (z != j)
						vectTemp.push_back(vectMatrix[i][z]);
				}
				vectMatrixS.push_back(vectTemp);
			}
			dSum += pow(-1, 0 - j) * vectMatrix[0][j] * MatrixDet(vectMatrixS, rows - 1);
		}
		return dSum;
	}
}

//矩阵求逆
vector> MatrixInv(const vector> &vectMatrix, int rows)
{
	vector> vectMatrixInv;
	vectMatrixInv = MatrixInit0(rows, rows);
	double dDet = MatrixDet(vectMatrix, rows);
	for (int i = 0; i < rows; i++) {
		for (int j = 0; j < rows; j++) {
			vector> vectTemp;
			vectTemp = MatrixInit0(rows - 1, rows - 1);
			int xi = 0, yj = 0;
			for (int g = 0; g < rows; g++) {
				for (int q = 0; q < rows; q++) {
					if (q != j && g != i) {
						vectTemp[xi][yj] = vectMatrix[g][q];
						yj++;
						if (yj == rows - 1) {
							xi++;
							yj = 0;
						}
					}
				}
			}
			double q;
			if (MatrixDet(vectTemp, rows - 1) == 0)
				q = 0;
			else
				q = pow(-1, i - j) * MatrixDet(vectTemp, rows - 1);
			vectMatrixInv[i][j] = q;
		}
	}
	vectMatrixInv = MatrixT(vectMatrixInv, rows, rows);
	for (int i = 0; i < rows; i++)
		for (int j = 0; j < rows; j++)
			vectMatrixInv[i][j] = vectMatrixInv[i][j] / dDet;

	return vectMatrixInv;
}

//矩阵相乘
vector> MatrixMul(const vector> &vectMatrixA, const vector> &vectMatrixB, 
	int rows1, int columns1, int rows2, int columns2)
{
	vector> vectMatrix;
	vectMatrix = MatrixInit0(rows1, columns2);
	for (int i = 0; i < rows1; i++) {
		for (int j = 0; j < columns2; j++) {
			double x = 0;
			for (int l = 0; l < columns1; l++)
				for (int r = 0; r < rows2; r++) {
					if (l == r)
						x += vectMatrixA[i][l] * vectMatrixB[r][j];
				}
			vectMatrix[i][j] = x;
		}
	}
	return vectMatrix;
}