矩陣類及其常規運算(加、減、乘、轉置、求逆、行列式、代數餘子式、伴隨矩陣)

PRE_MITER發表於2020-12-24

寫在前面

本程式中矩陣類的定義為師長給予(兩個SetData函式除外),在此不作詳細說明。其中涉及的對新舊資料的處理、對函式的封裝與保護、對指標的運用及其他理念,由讀者自行體會。此外,本程式中對於矩陣外框列印所採用的對齊方法有很大的侷限性,其餘演算法亦極為粗糙簡陋,若讀者能夠指出不當之處,當不勝感激。

程式碼實現

#include <iostream>
#include <iomanip>
#include <cstring>
#include <cmath>
using namespace std;

#define ull unsigned long long

class Matrix
{
	double* pData;
	int rowCount;
	int colCount;

public:
	Matrix(int rows, int cols) :
		pData(nullptr), rowCount(0), colCount(0)
	{
		CreateData(rows, cols);
		ClearData();
	}

	Matrix(int rows, int cols, double value) :
		pData(nullptr), rowCount(0), colCount(0)
	{
		CreateData(rows, cols);
		ClearData(value);
	}

	Matrix(int rows, int cols, const double* initArr, int initArrSize) :
		pData(nullptr), rowCount(0), colCount(0)
	{
		CreateData(rows, cols);
		CopyDataFrom(initArr, initArrSize);
	}

	~Matrix()
	{
		if (pData) delete[]pData;
	}

	Matrix(const Matrix& m) :
		pData(nullptr), rowCount(0), colCount(0)
	{
		CreateData(m.rowCount, m.colCount);
		CopyDataFrom(m.pData, m.rowCount * m.colCount);
	}

	Matrix& operator =(const Matrix& m)
	{
		CreateData(m.rowCount, m.colCount);
		CopyDataFrom(m.pData, m.rowCount * m.colCount);

		return *this;
	}

	double GetData(int rowIndex, int colIndex) const
	{
		return pData[rowIndex * colCount + colIndex];
	}

	double* SetData() const
	{
		return pData;
	}

private:
	void CreateData(int rows, int cols)
	{
		if (rows <= 0) rows = 1;
		if (cols <= 0) cols = 1;

		double* pOldData = nullptr;
		if (pData)
			if (rowCount * colCount != rows * cols)
				pOldData = pData;

		pData = new double[rows * cols];
		this->rowCount = rows;
		this->colCount = cols;
	}

	void CopyDataFrom(const double* newData, int newDataSize, bool clearOther = true)
	{
		double* pTarget = pData;
		const double* pSource = newData;

		double* pTargetEnd1 = pData + rowCount * colCount;
		double* pTargetEnd2 = pData + newDataSize;

		double* pTargetEnd = (pTargetEnd1 < pTargetEnd2) ? pTargetEnd1 : pTargetEnd2;

		while (pTarget < pTargetEnd)
			*(pTarget++) = *(pSource++);

		/*   for (; pTarget < pTargetEnd; ++pTarget, ++pSource)
			   *pTarget = *pSource;

		   for (int i = 0; i < newDataSize && i < rowCount * colCount; ++i)
			   pData[i] = newData[i];*/ //?上方迴圈的另外兩種寫法。

		if (clearOther)
			while (pTarget < pTargetEnd1)
				*pTarget++ = 0;
	}

	void ClearData(double value = 0)
	{
		double* pTarget = pData;
		double* pTargetEnd = pData + rowCount * colCount;

		while (pTarget < pTargetEnd)
			*(pTarget++) = value;
	}

	int GetRowCount() const { return rowCount; }
	int GetColCount() const { return colCount; }
	friend ostream& operator << (ostream& out, const Matrix& m);
	friend int GetRowCount(const Matrix& m);
	friend int GetColCount(const Matrix& m);
	friend double* SetData(Matrix& m, const int rowIndex, const int colIndex);
}; //矩陣類

int GetRowCount(const Matrix& m)
{
	return m.rowCount;
}

int GetColCount(const Matrix& m)
{
	return m.colCount;
}

//返回矩陣指向特定行列的指標
double* SetData(Matrix& m, const int rowIndex, const int colIndex) 
{
	double* p = m.pData;
	p += rowIndex * m.colCount + colIndex;
	return p;
}

//輸出流過載(使其能夠列印矩陣)?
ostream& operator << (ostream& out, const Matrix& m) 
{
	int x = 0, y = 0, i = 0; 
	ull blank = 0; ull numsize = 0;
	double temp = 0, edge = 0;

	blank = 9 * (ull)m.colCount - 1;
	out << "┏" << setw(blank) << "┓" << endl;
	for (y = 0; y < m.rowCount; ++y)
	{
		out << "┃";
		
		for (x = 0; x < m.colCount; ++x)
		{
			temp = m.GetData(y, x);  //儲存本次要列印的值
			if (temp < 1e-10) temp = 0; //temp過小直接設為0
			edge = fabs(temp - (long)temp); //判斷是否為整數

			if (temp != 0 && edge < 1e-4) edge = 0; //不為整數,將edge設為0,將其當作整數
			else if (edge >= 1e-4) edge = 1; //將其當作小數

			numsize = (ull)log10(fabs(temp)) + 1;
			if (temp < 0) numsize += 1;
			if (fabs(temp) < 1) numsize = 0;

			if (numsize == 0 && edge == 1) //處理小於1且帶小數的情況。
				out << setprecision(5) << setw(8) << temp;
			else if (numsize == 1)
			{
				if (edge == 1)
					out << setprecision(5) << setw(8) << temp;
				else
					out << setw(5) << temp << setw(3) << "";
			}
			else 
				out << setw(5) << temp << setw(3) << "";

			if (x == m.colCount - 1) out << "┃";
		}
		out << endl;
		out << (y == m.rowCount - 1 ? "┗" : "");
		if (y == m.rowCount - 1) out << setw(blank);
		out << (y == m.rowCount - 1 ? "┛" : "");}
	
	return out;
}

//矩陣轉置?
Matrix operator ~ (const Matrix& m) 
{
	Matrix mT(GetColCount(m), GetRowCount(m));
	double* p = mT.SetData();
	for (int row = 0; row < GetRowCount(mT); ++row)
	{
		for (int col = 0; col < GetColCount(mT); ++col)
			*p++ = m.GetData(col, row);
	}

	return mT;
}

//矩陣乘法?
Matrix operator * (const Matrix& m1, const Matrix& m2) 
{
	Matrix mX(GetRowCount(m1), GetColCount(m2));
	Matrix ErrormX(1, 1, 0);
	if (GetColCount(m1) != GetRowCount(m2))
	{
		cout << "兩矩陣行列不對應,無法相乘。返回錯誤值:1階0矩陣。" << endl;
		return ErrormX;
	}

	int t1 = GetColCount(m1), t2 = GetColCount(m2), t3 = GetRowCount(m1);
	int i = 0, x = 0, y = 0;

	double* pTarget = mX.SetData();

	while (x < t3) //t3為m1的行數
	{
		while (y < t2) //t2為m2的列數
		{
			while (i < t1) //t1為m1的列數(也是m2的行數),對每一個y(即m2的每一列),迴圈t1次計算出mX當前位置的值,然後y+1
			{
				*pTarget += m1.GetData(x, i) * m2.GetData(i, y);
				i++;
			}
			i = 0; y++; pTarget++;
		}
		y = 0; x++;
	}

	return mX;
}

//矩陣加法?
Matrix operator + (const Matrix& m1, const Matrix& m2)
{
	Matrix mAdd(GetRowCount(m1), GetColCount(m1));
	Matrix m_Error(1, 1, 0);
	double* p = mAdd.SetData();

	if (GetRowCount(m1) != GetRowCount(m2) ||
		GetColCount(m1) != GetColCount(m2))
	{
		cout << "兩矩陣行列不完全相同,無法相加!返回錯誤值:1階0矩陣。";
		return m_Error;
	}

	for (int row = 0; row < GetRowCount(m1); row++)
	{
		for (int col = 0; col < GetColCount(m1); col++)
			*p++ = m1.GetData(row, col) + m2.GetData(row, col);
	}

	return mAdd;
}

//矩陣減法
Matrix operator - (const Matrix& m1, const Matrix& m2)
{
	Matrix mMinus(GetRowCount(m1), GetColCount(m1));
	Matrix m_Error(1, 1, 0);
	double* p = mMinus.SetData();

	if (GetRowCount(m1) != GetRowCount(m2) ||
		GetColCount(m1) != GetColCount(m2))
	{
		cout << "兩矩陣行列不完全相同,無法相減!返回錯誤值:1階0矩陣。";
		return m_Error;
	}

	for (int row = 0; row < GetRowCount(m1); row++)
	{
		for (int col = 0; col < GetColCount(m1); col++)
			*p++ = m1.GetData(row, col) - m2.GetData(row, col);
	}

	return mMinus;
}


//交換函式,用於進行兩行交換
inline void swap(Matrix& mdet, int r1, int r2)
{
	int t = GetColCount(mdet);
	for (int i = 0; i < t; i++)
		swap(*SetData(mdet, r1, i), *SetData(mdet, r2, i));
}
//乘積函式,用於對矩陣指定行實行數乘
inline void Matrix_Multiline(Matrix& m, double k, int R)
{
	int C = GetColCount(m);
	for (int i = 0; i < C; i++)
		*SetData(m, R, i) *= k;
}

// 高斯消元求行列式
double fabs(Matrix& m)
{
	Matrix mDen = m;
	int row = GetRowCount(m), col = GetColCount(m);
	if (row != col) //判斷是否能求行列式。
	{
		cout << "該矩陣不是方陣,無法求行列式. 返回錯誤值0." << endl;
		cout << "◢";
		return 0;
	}

	if (row == 2)
		return m.GetData(0, 0) * m.GetData(1, 1) - m.GetData(0, 1) * m.GetData(1, 0);
	
	/*
	引數說明:
	1. k    : 儲存最終計算結果。
	2. temp : 當目前迴圈至t列時,儲存t+1行位於主對角線下方的數值。
	3. i,j  : 控制迴圈需要。
	4. count: 判斷首列是否全為0.
	5. t    : 控制迴圈沿對角線方向下降。
	*/
	double k = 1, temp = 0; int i = 0, j = 0, count = 0, t = 0;

	for (i = 0; i < row; i++)
	{
		if (mDen.GetData(i, 0) == 0)
			count += 1;
	}
	if (count == row) return 0;

	for (t = 0; t < col; t++)//從首列開始,將矩陣變為上三角。沿對角線方向下降,下文將該列主對角線上的元素稱為該列首項。
	{
		// 將該列首項變為1. ?
		if (mDen.GetData(t, t) == 1) {} //若該列首項為1,不做任何事。
		else
		{
			for (i = t; i < row; i++) //從該列首項開始。
			{
				if (mDen.GetData(i, t) == 1) //如果該列i行為1,將其與該列首項所在行對換,跳出行迴圈
				{
					if (i != t)
					{
						swap(mDen, i, t); //執行對換。
						k = (i == t ? k : -k); //行對換,符號改變。
					}
					goto Deal; //跳轉至處理函式。
				}
				else
				{
					if (mDen.GetData(i, t) == 0) //遇到0,進入下一行
						continue;
					else //如果該列i行不為0且不為1,用k儲存該行首項的值,並將該行首項變為1,之後將該行與該列首項所在行對換。
					{
						temp = mDen.GetData(i, t); //儲存該行首項的值。
						k *= temp; 
						//cout << "k = " << k << endl << endl;
						Matrix_Multiline(mDen, (1 / temp), i); //將該行首項變為1.
						if (i != t)
						{
							swap(mDen, i, t); //交換該行與該列首項所在行。
							k = -k; //改變符號。
						}
						goto Deal; //跳轉至處理函式。
					}
				}
			}
			if (t == 0 && count == row) return 0; // 首列全為0,返回0。
		}

	Deal: //將該列除首項外所有元素變為0.?
		for (i = t + 1; i < row; i++) //從該列首項所在下一行開始。
		{
			temp = mDen.GetData(i, t); //儲存該行首項的值。
			if (temp != 0) //如果不為0
			{
				for (j = t; j < col; j++) //從該行主對角線下項開始,每一項都減去將該列首項(即1)所在的行乘上temp。將該行首項變為0.
				{
					//cout << "第" << i + 1 << "行" << j + 1 << "列的元素為" << mDen.GetData(i, j) << endl; ?除錯用語句
					*SetData(mDen, i, j) -= temp * *SetData(mDen, t, j);
					//cout << "第" << i + 1 << "行" << j + 1 << "列的元素為" << mDen.GetData(i, j) << endl; ?除錯用語句
					//cout << "k = " << k << endl; ?除錯用語句
				}
			}
		}
	}

	return k;
}

// 提取餘子式對應矩陣
Matrix alogcom(Matrix& m, const int rowIndex, const int colIndex)
{
	int row = GetRowCount(m), col = GetColCount(m);
	int x = 0, y = 0, i = 0, j = 0;
	Matrix mC(row - 1, col - 1);

	for (y = 0; y < row; y++)
	{
		for (x = 0; x < col; x++)
		{
			if (y != rowIndex && x != colIndex)
			{
				*SetData(mC, i, j++) = *SetData(m, y, x);
				if (j == row - 1)
				{
					j = 0; i++;
					break;
				}
			}
		}
	}
	//cout << mC << endl << endl; ?除錯用語句
	return mC;
}
// 計算伴隨矩陣
Matrix adjoint(Matrix& m)
{
	int row = GetRowCount(m), col = GetColCount(m);
	short sign = -1;
	Matrix mCom(row, col);
	Matrix mTemp(row, col);
	for (int y = 0; y < row; y++)
	{
		for (int x = 0; x < col; x++)
		{
			mTemp = alogcom(m, x, y);
			//cout << "A" << x + 1 << y + 1 << " = " << fabs(mTemp) << endl; ?除錯用語句
			*SetData(mCom, y, x) = pow(sign, x + y + 2) * fabs(mTemp);
			//cout << mCom << endl; ?除錯用語句
		}
	}

	//cout << mCom << endl; ?除錯用語句
	return mCom;
}
// 求逆矩陣
Matrix operator !(Matrix& m)
{
	int row = GetRowCount(m), col = GetColCount(m);
	Matrix m_(row, col);
	Matrix mCom = adjoint(m);
	Matrix Errorm(1, 1, 0);
	double Den = fabs(m);
	//cout << Den << endl; ?除錯用語句
	if (row != col || Den == 0)
	{
		cout << "該矩陣不可逆。返回錯誤值:1階0矩陣。" << endl;
		return Errorm;
	}

	for (int y = 0; y < row; y++)
	{
		for (int x = 0; x < col; x++)
		{
			*SetData(m_, y, x) = mCom.GetData(y, x) / Den;
			//cout << m_ << endl; ?除錯用語句
		}
	}


	return m_;
}

int main()
{
	double arr[15] = { 1, 2, 3, 4, 5, 6, 1, 3, 4, 2, 3, 1, 5, 0.1, -0.5 };
	double arr_[9] = { 1, 0, 0, 0, 1, 0, 0, 0, 1 };
	Matrix m1(3, 3, arr, 9);
	Matrix m2 = ~m1;
	Matrix m3 = m1 * ~m1;
	Matrix m4(3, 3, arr_, 9);
	Matrix m5 = !m1;
	Matrix m6 = m1 * !m1;
	Matrix m7 = adjoint(m1);

	cout << "m1 = " << endl << m1 << endl << endl;
	cout << "m1^T(轉置) = " << endl << m2 << endl << endl;
	cout << "m1 * m1^T = " << endl << m3 << endl << endl;
	cout << "m1^*(伴隨) = " << endl << m7 << endl << endl;
	cout << "m1 * m1^* = "  << endl << m1 * m7 << endl << endl;
	cout << "|m1|(行列式) = " << fabs(m1) << endl << endl;
	cout << "m1^-1(逆矩陣) = " << endl << m5 << endl << endl;
	cout << "m1 * m1^-1 = " << endl << m6 << endl << endl;
	cout << "m4 = " << endl << m4 << endl << endl;
	cout << "|m4| = " << fabs(m4) << endl << endl;
	cout << "m1 + m4 = " << endl << m1 + m4 << endl << endl;
	cout << "m1 - m4 = " << endl << m1 - m4 << endl << endl;

	system("pause");
	return 0;
}

執行結果

在這裡插入圖片描述

寫在後面

矩陣的其他運算,等高等代數學明白後會補上?

相關文章