天天看點

在C++中實作矩陣運算

在C++中實作矩陣運算

本文部落格連結:http://blog.csdn.net/jdh99,作者:jdh,轉載請注明.

參考連結:http://sc.dhu.edu.cn/weblearning/math/jisuanfangfa/matlabpg/pg6.htm

環境:

主機:XP

開發環境:mingw

功能:

在C++中實作矩陣運算

說明:

定義了兩個類,一個是_Matrix,這是一個二維矩陣類,定義了初始化的方法.另一個是_Matrix_Calc,這個類定義了矩陣的基本運算,包括加,減,乘,轉置,求行列式值,求逆矩陣.

源代碼:

_Matrix.h:

#ifndef _MATRIX_H
#define _MATRIX_H

//頭檔案
#include <stdio.h>
#include <stdlib.h>

//矩陣資料結構
//二維矩陣
class _Matrix
{
private:
public:
	int m;
	int n;
	float *arr;
	
	//初始化
	_Matrix(int mm = 0,int nn = 0);
	//設定m
	void set_m(int mm);
	//設定n
	void set_n(int nn);
	//初始化
	void init_matrix();
	//釋放
	void free_matrix();
	//讀取i,j坐标的資料
	//失敗傳回-31415,成功傳回值
	float read(int i,int j);
	//寫入i,j坐标的資料
	//失敗傳回-1,成功傳回1
	int write(int i,int j,float val);
};

//二維運算類
class _Matrix_Calc
{
private:
public:
	_Matrix_Calc();
	//C = A + B
	//成功傳回1,失敗傳回-1
	int add(_Matrix *A,_Matrix *B,_Matrix *C);
	//C = A - B
	//成功傳回1,失敗傳回-1
	int subtract(_Matrix *A,_Matrix *B,_Matrix *C);
	//C = A * B
	//成功傳回1,失敗傳回-1
	int multiply(_Matrix *A,_Matrix *B,_Matrix *C);
	//行列式的值,隻能計算2 * 2,3 * 3
	//失敗傳回-31415,成功傳回值
	float _Matrix_Calc::det(_Matrix *A);
	//求轉置矩陣,B = AT
	//成功傳回1,失敗傳回-1
	int transpos(_Matrix *A,_Matrix *B);
	//求逆矩陣,B = A^(-1)
	//成功傳回1,失敗傳回-1
	int inverse(_Matrix *A,_Matrix *B);
};

#endif
           

_Matrix.cpp:

#include "_Matrix.h"

//矩陣類方法
//初始化
_Matrix::_Matrix(int mm,int nn)
{
	m = mm;
	n = nn;
}

//設定m
void _Matrix::set_m(int mm)
{
	m = mm;
}

//設定n
void _Matrix::set_n(int nn)
{
	n = nn;
}

//初始化
void _Matrix::init_matrix()
{
	arr = new float[m * n];
}

//釋放
void _Matrix::free_matrix()
{
	delete []arr;
}

//讀取i,j坐标的資料
//失敗傳回-31415,成功傳回值
float _Matrix::read(int i,int j)
{
	if (i >= m || j >= n)
	{
		return -31415;
	}
	
	return *(arr + i * n + j);
}

//寫入i,j坐标的資料
//失敗傳回-1,成功傳回1
int _Matrix::write(int i,int j,float val)
{
	if (i >= m || j >= n)
	{
		return -1;
	}
	
	*(arr + i * n + j) = val;
	return 1;
}

//矩陣運算類方法
//初始化
_Matrix_Calc::_Matrix_Calc()
{
}

//C = A + B
//成功傳回1,失敗傳回-1
int _Matrix_Calc::add(_Matrix *A,_Matrix *B,_Matrix *C)
{
	int i = 0;
	int j = 0;
	
	//判斷是否可以運算
	if (A->m != B->m || A->n != B->n || \
		A->m != C->m || A->n != C->n)
	{
		return -1;
	}
	//運算
	for (i = 0;i < C->m;i++)
	{
		for (j = 0;j < C->n;j++)
		{
			C->write(i,j,A->read(i,j) + B->read(i,j));
		}
	}
	
	return 1;
}

//C = A - B
//成功傳回1,失敗傳回-1
int _Matrix_Calc::subtract(_Matrix *A,_Matrix *B,_Matrix *C)
{
	int i = 0;
	int j = 0;
	
	//判斷是否可以運算
	if (A->m != B->m || A->n != B->n || \
		A->m != C->m || A->n != C->n)
	{
		return -1;
	}
	//運算
	for (i = 0;i < C->m;i++)
	{
		for (j = 0;j < C->n;j++)
		{
			C->write(i,j,A->read(i,j) - B->read(i,j));
		}
	}
	
	return 1;
}

//C = A * B
//成功傳回1,失敗傳回-1
int _Matrix_Calc::multiply(_Matrix *A,_Matrix *B,_Matrix *C)
{
	int i = 0;
	int j = 0;
	int k = 0;
	float temp = 0;
	
	//判斷是否可以運算
	if (A->m != C->m || B->n != C->n || \
		A->n != B->m)
	{
		return -1;
	}
	//運算
	for (i = 0;i < C->m;i++)
	{
		for (j = 0;j < C->n;j++)
		{
			temp = 0;
			for (k = 0;k < A->n;k++)
			{
				temp += A->read(i,k) * B->read(k,j);
			}
			C->write(i,j,temp);
		}
	}
	
	return 1;
}

//行列式的值,隻能計算2 * 2,3 * 3
//失敗傳回-31415,成功傳回值
float _Matrix_Calc::det(_Matrix *A)
{
	float value = 0;
	
	//判斷是否可以運算
	if (A->m != A->n || (A->m != 2 && A->m != 3))
	{
		return -31415;
	}
	//運算
	if (A->m == 2)
	{
		value = A->read(0,0) * A->read(1,1) - A->read(0,1) * A->read(1,0);
	}
	else
	{
		value = A->read(0,0) * A->read(1,1) * A->read(2,2) + \
				A->read(0,1) * A->read(1,2) * A->read(2,0) + \
				A->read(0,2) * A->read(1,0) * A->read(2,1) - \
				A->read(0,0) * A->read(1,2) * A->read(2,1) - \
				A->read(0,1) * A->read(1,0) * A->read(2,2) - \
				A->read(0,2) * A->read(1,1) * A->read(2,0);
	}
	
	return value;
}

//求轉置矩陣,B = AT
//成功傳回1,失敗傳回-1
int _Matrix_Calc::transpos(_Matrix *A,_Matrix *B)
{
	int i = 0;
	int j = 0;
	
	//判斷是否可以運算
	if (A->m != B->n || A->n != B->m)
	{
		return -1;
	}
	//運算
	for (i = 0;i < B->m;i++)
	{
		for (j = 0;j < B->n;j++)
		{
			B->write(i,j,A->read(j,i));
		}
	}
	
	return 1;
}

//列印2維矩陣
void printff_matrix(_Matrix *A)
{
	int i = 0;
	int j = 0;
	int m = 0;
	int n = 0;
	
	m = A->m;
	n = A->n;
	for (i = 0;i < m;i++)
	{
		for (j = 0;j < n;j++)
		{
			printf("%f ",A->read(i,j));
		}
		printf("\n");
	}
}

//求逆矩陣,B = A^(-1)
//成功傳回1,失敗傳回-1
int _Matrix_Calc::inverse(_Matrix *A,_Matrix *B)
{
	int i = 0;
	int j = 0;
	int k = 0;
	_Matrix m(A->m,2 * A->m);
	float temp = 0;
	float b = 0;
	
	//判斷是否可以運算
	if (A->m != A->n || B->m != B->n || A->m != B->m)
	{
		return -1;
	}
	
	/*
	//如果是2維或者3維求行列式判斷是否可逆
	if (A->m == 2 || A->m == 3)
	{
		if (det(A) == 0)
		{
			return -1;
		}
	}
	*/
	
	//增廣矩陣m = A | B初始化
	m.init_matrix();
	for (i = 0;i < m.m;i++)
	{
		for (j = 0;j < m.n;j++)
		{
			if (j <= A->n - 1)
			{
				m.write(i,j,A->read(i,j));
			}
			else
			{
				if (i == j - A->n)
				{
					m.write(i,j,1);
				}
				else
				{
					m.write(i,j,0);
				}
			}
		}
	}
	
	//高斯消元
	//變換下三角
	for (k = 0;k < m.m - 1;k++)
	{
		//如果坐标為k,k的數為0,則行變換
		if (m.read(k,k) == 0)
		{
			for (i = k + 1;i < m.m;i++)
			{
				if (m.read(i,k) != 0)
				{
					break;
				}
			}
			if (i >= m.m)
			{
				return -1;
			}
			else
			{
				//交換行
				for (j = 0;j < m.n;j++)
				{
					temp = m.read(k,j);
					m.write(k,j,m.read(k + 1,j));
					m.write(k + 1,j,temp);
				}
			}
		}
		
		//消元
		for (i = k + 1;i < m.m;i++)
		{
			//獲得倍數
			b = m.read(i,k) / m.read(k,k);
			//行變換
			for (j = 0;j < m.n;j++)
			{
				temp = m.read(i,j) - b * m.read(k,j);
				m.write(i,j,temp);
			}
		}
	}
	//變換上三角
	for (k = m.m - 1;k > 0;k--)
	{
		//如果坐标為k,k的數為0,則行變換
		if (m.read(k,k) == 0)
		{
			for (i = k + 1;i < m.m;i++)
			{
				if (m.read(i,k) != 0)
				{
					break;
				}
			}
			if (i >= m.m)
			{
				return -1;
			}
			else
			{
				//交換行
				for (j = 0;j < m.n;j++)
				{
					temp = m.read(k,j);
					m.write(k,j,m.read(k + 1,j));
					m.write(k + 1,j,temp);
				}
			}
		}
		
		//消元
		for (i = k - 1;i >= 0;i--)
		{
			//獲得倍數
			b = m.read(i,k) / m.read(k,k);
			//行變換
			for (j = 0;j < m.n;j++)
			{
				temp = m.read(i,j) - b * m.read(k,j);
				m.write(i,j,temp);
			}
		}
	}
	//将左邊方陣化為機關矩陣
	for (i = 0;i < m.m;i++)
	{
		if (m.read(i,i) != 1)
		{
			//獲得倍數
			b = 1 / m.read(i,i);
			//行變換
			for (j = 0;j < m.n;j++)
			{
				temp = m.read(i,j) * b;
				m.write(i,j,temp);
			}
		}
	}
	//求得逆矩陣
	for (i = 0;i < B->m;i++)
	{
		for (j = 0;j < B->m;j++)
		{
			B->write(i,j,m.read(i,j + m.m));
		}
	}
	//釋放增廣矩陣
	m.free_matrix();
	
	return 1;
}
           

main.cpp:(測試代碼)

#include <stdio.h>
#include <iostream>
#include "_Matrix.h"

//帶速度變量卡爾曼濾波
using namespace std;

//列印2維矩陣
void printf_matrix(_Matrix *A)
{
	int i = 0;
	int j = 0;
	int m = 0;
	int n = 0;
	
	m = A->m;
	n = A->n;
	for (i = 0;i < m;i++)
	{
		for (j = 0;j < n;j++)
		{
			printf("%f\t",A->read(i,j));
		}
		printf("\n");
	}
}

int main()
{
	int i = 0;
	int j = 0;
	int k = 0;
	_Matrix_Calc m_c;
	_Matrix m1(3,3);
	_Matrix m2(3,3);
	_Matrix m3(3,3);
	
	//初始化記憶體
	m1.init_matrix();
	m2.init_matrix();
	m3.init_matrix();
	
	//初始化資料
	k = 1;
	for (i = 0;i < m1.m;i++)
	{
		for (j = 0;j < m1.n;j++)
		{
			m1.write(i,j,k++);
		}
	}
	
	for (i = 0;i < m2.m;i++)
	{
		for (j = 0;j < m2.n;j++)
		{
			m2.write(i,j,k++);
		}
	}
	
	//原資料
	printf("A:\n");
	printf_matrix(&m1);
	printf("B:\n");
	printf_matrix(&m2);
	
	printf("A:行列式的值%f\n",m_c.det(&m1));
	
	//C = A + B
	if (m_c.add(&m1,&m2,&m3) > 0)
	{
		printf("C = A + B:\n");
		printf_matrix(&m3);
	}
	
	//C = A - B
	if (m_c.subtract(&m1,&m2,&m3) > 0)
	{
		printf("C = A - B:\n");
		printf_matrix(&m3);
	}
	
	//C = A * B
	if (m_c.multiply(&m1,&m2,&m3) > 0)
	{
		printf("C = A * B:\n");
		printf_matrix(&m3);
	}
	
	//C = AT
	if (m_c.transpos(&m1,&m3) > 0)
	{
		printf("C = AT:\n");
		printf_matrix(&m3);
	}
	
	/*
	m1.write(0,0,0);
	m1.write(0,1,1);
	m1.write(0,2,0);
	m1.write(1,0,1);
	m1.write(1,1,1);
	m1.write(1,2,-1);
	m1.write(2,0,0);
	m1.write(2,1,-2);
	m1.write(2,2,3);
	*/
	if (m_c.inverse(&m1,&m3) > 0)
	{
		printf("C = A^(-1):\n");
		printf_matrix(&m3);
	}
	
	getchar();
	return 0;
}