原理基於 gaussian jordan elimination 方法，考慮求解如下線性方程組：

$$

\begin{equations}

\mathbf{A}\mathbf{x} = \mathbf{b}

\end{equations}

$$

　　我們設 $ L $、$ U $、$ P $ 滿足：

$$

\begin{equations}
\mathbf{PA}&=\mathbf{LU}
\end{equations}

$$

　　所以只要找到這樣的$ L $、$ U $、$ P $， 就可以通過：

$$

\begin{equations}

\mathbf{Ax}&=\mathbf{b} \\
\mathbf{PAx}&=\mathbf{Pb} \\
\mathbf{LUx}&=\mathbf{Pb} \\
\mathbf{Ux}&=\mathbf{y} \\
\mathbf{Ly}&=\mathbf{Pb}

\end{equations}

$$

　　求出 $ x $ 了，注意這裡 $ x $ 是一個列向量。所以該方法可以用來求矩陣的逆。

　　這裡 $ L $、$ U $、$ P $ 都是我們在做 gaussian jordan elimination 方法時交換行所產生的矩陣（程式碼中 ```lup_decomposition``` 函式中實現了該過程。）

　　程式碼實現：

#include <vector>
#include <cassert>
#include <ostream>
#include <iomanip>
#include <iostream>

using std::vector;
using std::ostream;
using std::cout;
using std::swap;

bool lup_decomposition(vector<vector<double>> a, vector<vector<double>>& l, vector<vector<double>>& u, vector<double>& p)
{
	for (int i = 0; i < a.size(); i++)
		p[i] = (double)i;
	for (int i = 0; i < a.size(); i++)
	{
		l[i][i] = 1;
		for (int j = i + 1; j < a.size(); j++)
			l[i][j] = 0;
	}
	for (int i = 1; i < a.size(); i++)
		for (int j = 0; j < i; j++)
			u[i][j] = 0;
	for (int i = 0; i < a.size(); i++)
	{
		double n = 0; int i_;
		for (int j = i; j < a.size(); j++)
			if (abs(a[j][i]) > n)
				n = abs(a[j][i]), i_ = j;
		if (n == 0) return false;
		swap(p[i], p[i_]);
		for (int j = 0; j < a.size(); j++)
			swap(a[i][j], a[i_][j]);
		u[i][i] = a[i][i];
		for (int j = i + 1; j < a.size(); j++)
		{
			l[j][i] = a[j][i] / u[i][i];
			u[i][j] = a[i][j];
		}
		for (int j = i + 1; j < a.size(); j++)
			for (int k = i + 1; k < a.size(); k++)
				a[j][k] -= l[j][i] * u[i][k];
	}
	return true;
}

void lup_solve(vector<vector<double>> a, vector<double> b, vector<vector<double>> l, vector<vector<double>> u, vector<double> p, vector<double>& x)
{
	vector<double> y(a.size(), 0);
	for (int i = 0; i < a.size(); i++)
	{
		double sum = 0;
		for (int j = 0; j < i; j++)
			sum += l[i][j] * y[j];
		y[i] = b[p[i]] - sum;
	}
	for (int i = a.size() - 1; i >= 0; i--)
	{
		double sum = 0;
		for (int j = i + 1; j < a.size(); j++)
			sum += u[i][j] * x[j];
		x[i] = (y[i] - sum) / u[i][i];
	}
}

vector<vector<double>> inverse(vector<vector<double>> a)
{
	vector<vector<double>>
		unit_matrix(a.size(), vector<double>(a.size(), 0)),
		x(a.size(), vector<double>(a.size(), 0)),
		l(a.size(), vector<double>(a.size(), 0)),
		u(a.size(), vector<double>(a.size(), 0));
	vector<double> p(a.size());
	for (int i = 0; i < a.size(); i++)
		unit_matrix[i][i] = 1;
	lup_decomposition(a, l, u, p);
	for (int i = 0; i < a.size(); i++)
		lup_solve(a, unit_matrix[i], l, u, p, x[i]);
	return x;
}

vector<double> sum(vector<vector<double>> a, int axis = 1)
{
	if (axis == 1)
	{
		vector<double> e(a.size());
		for (int i = 0; i < a.size(); i++)
			for (int j = 0; j < a[0].size(); j++)
				e[i] += a[i][j];
		return e;
	}
	else if (axis == -1)
	{
		vector<double> e(a[0].size());
		for (int i = 0; i < a[0].size(); i++)
			for (int j = 0; j < a.size(); j++)
				e[i] += a[i][j];
		return e;
	}
	else {
		vector<double> e(1);
		for (int i = 0; i < a.size(); i++)
			for (int j = 0; j < a[0].size(); j++)
				e[0] += a[i][j];
		return e;
	}
}

double sum(vector<double> a)
{
	double e = 0;
	for (auto i : a)
		e += i;
	return e;
}

vector<double> add(vector<double> a, vector<double> b)
{
	assert(a.size() == b.size());
	vector<double> c(a.size());
	for (int i = 0; i < a.size(); i++)
		c[i] = a[i] + b[i];
	return c;
}

vector<vector<double>> add(vector<vector<double>> a, vector<vector<double>> b)
{
	assert(a.size() == b.size());
	assert(a[0].size() == b[0].size());
	vector<vector<double>> c(a.size(), vector<double>(a[0].size(), 0));
	for (int i = 0; i < a.size(); i++)
		for (int j = 0; j < a[0].size(); j++)
			c[i][j] = a[i][j] + b[i][j];
	return c;
}

vector<double> sub(vector<double> a, vector<double> b)
{
	assert(a.size() == b.size());
	vector<double> c(a.size());
	for (int i = 0; i < a.size(); i++)
		c[i] = a[i] - b[i];
	return c;
}

vector<vector<double>> sub(vector<vector<double>> a, vector<vector<double>> b)
{
	assert(a.size() == b.size());
	assert(a[0].size() == b[0].size());
	vector<vector<double>> c(a.size(), vector<double>(a[0].size(), 0));
	for (int i = 0; i < a.size(); i++)
		for (int j = 0; j < a[0].size(); j++)
			c[i][j] = a[i][j] - b[i][j];
	return c;
}

vector<double> mul(vector<vector<double>> a, vector<double> b)
{
	assert(a[0].size() == b.size());
	vector<double> c(b.size(), 0);
	for (int i = 0; i < a.size(); i++)
		for (int j = 0; j < a[0].size(); j++)
			c[i] += a[i][j] * b[j];
	return c;
}

vector<vector<double>> mul(vector<vector<double>> a, vector<vector<double>> b)
{
	assert(a[0].size() == b.size());
	vector<vector<double>> r(a.size(), vector<double>(b[0].size(), 0));
	for (int i = 0; i < r.size(); i++)
		for (int j = 0; j < r[0].size(); j++)
			for (int k = 0; k < b.size(); k++)
				r[i][j] += a[i][k] * b[k][j];
	return r;
}

double dot(vector<double> a, vector<double> b)
{
	assert(a.size() == b.size());
	double c = 0.;
	for (int i = 0; i < a.size(); i++)
		c += a[i] * b[i];
	return c;
}

vector<vector<double>> transpose(vector<vector<double>> a)
{
	for (int i = 0; i < a.size(); i++)
		for (int j = i + 1; j < a[0].size(); j++)
			swap(a[i][j], a[j][i]);
	return a;
}

ostream& operator<<(ostream& os, vector<vector<double>> a)
{
	os << "[\n";
	for (auto rows : a)
	{
		for (auto col : rows)
			os << std::right << std::setw(12) << col << ' ';
		os << '\n';
	}
	os << "]\n";
	return os;
}

ostream& operator<<(ostream& os, vector<double> a)
{
	for (auto i : a) os << i << '\n';
	return os;
}

int main()
{
	vector<vector<double>> m{ {1, 2}, {2, 3} };

	cout << m << '\n';
	cout << transpose(inverse(m)) << '\n';
	
	return 0;
}

　　測試結果：

　　當然程式碼也還有可以優化的地方，以及這些都是爭對方陣而編寫的一段矩陣、向量運算相關的簡陋程式碼。