Код для обращения таков (взят
отсюда):
Код:
// Обращение квадратной матрицы
template<typename TT>
inline static bool InvertMN(int N, TT **M,
int *GausOstatok=NULL,
long double *GausDet = NULL,
long double *GausMinved= NULL) {
// (c) Drobotenko http://www.drobotenko.com/bib_rus.html
int gaus_ostatok; // дефект линейного преобразования 0 == ОК
long double gaus_deter; // определитель в случае успешного обращения
long double gaus_minved; // минимальный ведущий элемент
// можно использовать для оценки точности
int *rn, *cn;
int j, k;
#define _NaN() (double&)*"Implementation dependent";
// заполнитель для неразрешенных измерений
// при невозможности обращения
rn = new int [N];
cn = new int [N];
for(j=N; j--; )
rn[j]=cn[j]=j;
gaus_minved=1e99;
gaus_deter=1;
for(gaus_ostatok=N; gaus_ostatok; gaus_ostatok--) {
int jved,kved;
long double vved=-1, t;
for(j=N; j--; ) { // поиск ведущего
if(~rn[j])
for(k=N; k--; )
if(~cn[k])
if(vved<fabsl(M[j][k]))
vved=fabsl(M[j][k]),jved=j,kved=k; }
if(gaus_minved>vved)
gaus_minved=vved;
gaus_deter*=M[jved][kved];
if(vved<1e-99) { // Для long double можно уменьшить
for(j=N; j--; ) {
if(~rn[j])
for(k=N; k--; )
M[j][k]=_NaN(); // ?
if(~cn[j])
for(k=N; k--; )
M[k][j]=_NaN(); }
delete[] rn;
delete[] cn;
return false; }
int jt=rn[jved],kt=cn[kved];
// перестановки
for(j=N; j--; )
t=M[kt][j],M[kt][j]=M[jved][j],M[jved][j]=t;
for(j=N; j--; )
t=M[j][jt],M[j][jt]=M[j][kved],M[j][kved]=t;
rn[jved]=rn[kt];
cn[kved]=cn[jt];
rn[kt]=cn[jt]=-1;
vved=M[kt][jt]; M[kt][jt]=1;
for(j=N; j--; ) {
if(j==kt)
continue;
long double mul=M[j][jt]/vved;
M[j][jt]=0;
for(k=N; k--; )
M[j][k]-=M[kt][k]*mul; } // самый внутренний цикл ровно N^3 операций
for(k=N; k--; )
M[kt][k]/=vved;
}
if (GausOstatok) *GausOstatok=gaus_ostatok;
if (GausDet) *GausDet = gaus_deter;
if (GausMinved) *GausMinved= gaus_minved;
delete[] rn;
delete[] cn;
return true;
}
Код работоспособный, при перемножении матриц

размером 3x3 погрешности имеют порядок -19.
Матрицы, используемые в выражениях, заданы, предположительно, верно.
Выражение

в случае

вычисляется так:

,
аналогичные вычисления в MatLab дают те же результаты.