Odwracanie macierzy to poważny problem (szczególnie dla studentów). Z tego powodu udostępniamy nasz kod C++ pozwalający na łatwe, szybkie i przyjemne odwrócenie macierzy kwadratowej o różnym wymiarze. Zaimplementowane zostały dwa popularne algorytmy: Gaussian elimination oraz LU decomposition(factorization) . Z oczywistych względów polecamy korzystanie z tego drugiego. Użycie jest bardzo łatwe. Wystarczy dowolną tablicę o rozmiarze N * N (gdzie N jest rozmiarem wiersza lub kolumny) , w postaci wskaźnika przekazać w sposób następujący:
<typ> X[N * N]; // gdzie <typ> to double/float/int/...
sophisticatedAlgebra::Inverter<typ, N>::gaussian(X); //LU
//albo
sophisticatedAlgebra::Inverter<typ, N>::lu(X); //LU
Obliczenia zostaną wykonane w miejscu, więc zawartość tablicy X zostanie nadpisana.
Oto nagłówek, zawierający niezbędny kod:
#ifndef INVERTER_H
#define INVERTER_H
#include <cmath>
namespace sophisticatedAlgebra
{
template <typename T>
T sgn(const T& val)
{
if (val >= static_cast<T>(0))
{
return static_cast<T>(1);
}
else
{
return static_cast<T>(-1);
}
}
template <typename T, int N>
class Inverter
{
public:
static void gaussian(T* X)
{
T a[(2 * N) + 1][(2 * N) + 1];
T d = 0;
for (int i = 0; i < 2 * N + 1; ++i)
{
for (int j = 0; j < 2 * N + 1; ++j)
{
a[i][j] = 0;
}
}
for (int i = 1; i <= N; i++)
{
for (int j = 1; j <= 2 * N; j++)
{
if (j == (i + N))
{
a[i][j] = 1;
}
}
}
for (int i = 1; i <= N; ++i)
{
for (int j = 1; j <= N; ++j)
{
a[i][j] = X[j - 1 + (i - 1) * N];
}
}
for (int i = N; i > 1; --i)
{
if (a[i - 1][1] < a[i][1])
{
for (int j = 1; j <= N * 2; ++j)
{
d = a[i][j];
a[i][j] = a[i - 1][j];
a[i - 1][j] = d;
}
}
}
for (int i = 1; i <= N; ++i)
{
for (int j = 1; j <= N * 2; ++j)
{
if (j != i)
{
d = a[j][i] / a[i][i];
for (int k = 1; k <= N * 2; ++k)
{
a[j][k] = a[j][k] - (a[i][k] * d);
}
}
}
}
for (int i = 1; i <= N; ++i)
{
d = a[i][i];
for (int j = 1; j <= N * 2; ++j)
{
a[i][j] /= d;
}
}
for (int i = 1; i <= N; ++i)
{
for (int j = N + 1; j <= N * 2; ++j)
{
X[j - (N + 1) + (i - 1) * N] = a[i][j];
}
}
}
static void lu(T* X)
{
T Z[N][N];
T lower[N][N];
T upper[N][N];
T I[N][N];
for (int column = 0; column < N; ++column)
{
for (int row = 0; row < N; ++row)
{
Z[column][row] = 0;
if (column == row)
{
I[column][row] = 1;
}
else
{
I[column][row] = 0;
}
}
}
luDecomposition(X, lower, upper);
for (int column = 0; column < N; ++column)
{
for (int row = 0; row < N; ++row)
{
Z[row][column] = valZ(lower, upper, I, Z, row, column);
}
}
for (int column = 0; column < N; ++column)
{
for (int row = N - 1; row >= 0; --row)
{
X[column + row * N] = inversal(upper, X, Z, row, column);
}
}
}
private:
static void luDecomposition(T *X, T lower[N][N], T upper[N][N])
{
for (int row = 0; row < N; ++row)
{
for (int column = 0; column < N; ++column)
{
if (column < row)
{
lower[column][row] = 0;
}
else
{
lower[column][row] = X[row + column * N];
for (int k = 0; k < row; ++k)
{
lower[column][row] = lower[column][row] - lower[column][k] * upper[k][row];
}
}
}
for (int column = 0; column < N; ++column)
{
if (column < row)
{
upper[row][column] = 0;
}
else if (column == row)
{
upper[row][column] = 1;
}
else
{
upper[row][column] = X[row + column * N] / lower[row][row];
for (int k = 0; k < row; k++)
{
upper[row][column] = upper[row][column] -
((lower[row][k] * upper[k][column]) / lower[row][row]);
}
}
}
}
}
static T valZ(const T lower[N][N],
const T upper[N][N],
const T I[N][N],
const T Z[N][N],
int row, int column)
{
T sum = 0;
for (int i = 0; i < N; i++)
{
if (i != row) {
sum += lower[row][i] * Z[i][column];
}
}
T result = I[row][column] - sum;
result = result / lower[row][row];
return result;
}
static T inversal(const T upper[N][N],
const T* inversed,
const T Z[N][N],
int row, int column)
{
T sum = 0;
for (int i = 0; i < N; i++)
{
if (i != row)
{
sum += upper[row][i] * inversed[column + i * N];
}
}
T result = Z[row][column] - sum;
result = result / upper[row][row];
return result;
}
};
template <typename T>
class Inverter<T, 0>
{
static void gaussian(T* X) {}
static void lu(T* X) {}
};
template <typename T>
class Inverter<T, 1>
{
static void gaussian(T* X)
{
X[0] = static_cast<T>(1) / X[0];
}
static void lu(T* X)
{
gaussian(X);
}
};
template <typename T>
class Inverter<T, 2>
{
static void gaussian(T* X)
{
T determinant;
T adjoint[4];
adjoint[0] = X[3];
adjoint[1] = -X[1];
adjoint[2] = -X[2];
adjoint[3] = X[0];
determinant = X[0] * adjoint[0] + X[1] * adjoint[2];
determinant = static_cast<T>(1) / determinant;
X[0] = adjoint[0] * determinant;
X[1] = adjoint[1] * determinant;
X[2] = adjoint[2] * determinant;
X[3] = adjoint[3] * determinant;
}
static void lu(T* X)
{
gaussian(X);
}
};
template <typename T>
class Inverter<T, 3>
{
static void gaussian(T* X)
{
T determinant;
T adjoint[9];
adjoint[0] = X[4] * X[8] - X[5] * X[7];
adjoint[1] = -X[1] * X[8] + X[2] * X[7];
adjoint[2] = X[1] * X[5] - X[2] * X[4];
adjoint[3] = -X[3] * X[8] + X[5] * X[6];
adjoint[4] = X[0] * X[8] - X[2] * X[6];
adjoint[5] = -X[0] * X[5] + X[2] * X[3];
adjoint[6] = X[3] * X[7] - X[4] * X[6];
adjoint[7] = -X[0] * X[7] + X[1] * X[6];
adjoint[8] = X[0] * X[4] - X[1] * X[3];
determinant = X[0] * adjoint[0] + X[1] * adjoint[3] + X[2] * adjoint[6];
determinant = static_cast<T>(1) / determinant;
X[0] = adjoint[0] * determinant;
X[1] = adjoint[1] * determinant;
X[2] = adjoint[2] * determinant;
X[3] = adjoint[3] * determinant;
X[4] = adjoint[4] * determinant;
X[5] = adjoint[5] * determinant;
X[6] = adjoint[6] * determinant;
X[7] = adjoint[7] * determinant;
X[8] = adjoint[8] * determinant;
}
static void lu(T* X)
{
gaussian(X);
}
};
}
#endif
Miłej laborki 
PS. Jeżeli wersja Gaussian elimination produkuje błędne wyniki, to jest to zwykle błąd numeryczny, na który ta metoda jest bardzo wrażliwa. To jeden z wielu powodów żeby wybrać LU decomposition.