Привет! Интересует метод нахождения обратной матрицы методом вращений с QR-разложением. Если у кого есть такой плз помогите. Ниже мой вариант. Еще очень нужен этот же метод, но блочный. Т.е. разбиение матрицы на блоки матрицы, для того чтобы каждая матрица помешалась в кеш, для ускорения работы, естественно с использованием данного метода.
#include <math.h>
#include "inv.h"
#define Eps 1.0e-100
int InvMatrix(int n, double *a, double *x)
{
int i, j, k;
double tmp, q, t, x1, x2;
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
x[i*n+j] = (i == j);

for (i = 0; i < n; i++)
{
for (j = n — 1; j > i; j--)
{
if (fabs(a[j*n+i]) < Eps)
continue;

tmp = a[(j-1)*n+i];
q = a[n*j+i];
t = sqrt(tmp * tmp + q * q);
tmp = tmp/t;
q = -q/t;

for (k = i; k < n; k++)
{
x1 = a[(j — 1)*n+k];
x2 = a[j*n+k];
a[(j — 1)*n+k] = x1 * tmp — x2 * q;
a[j*n+k] = x1 * q + x2 * tmp;
}
for (k = 0; k < n; k++)
{
x1 = x[k*n+(j — 1)];
x2 = x[k*n+j];
x[k*n+(j — 1)] = x1 * tmp — x2 * q;
x[k*n+j] = x1 * q + x2 * tmp;
}
}

if (fabs(a[i*n+i]) < Eps)
return 0;
}
for (i = 0; i < n; i++)
for (j = i + 1; j < n; j++)
{
tmp = x[i*n+j];
x[i*n+j] = x[j*n+ i];
x[j*n+i] = tmp;
}
for (k = 0; k < n; k++)
for (i = n — 1; i >= 0; i--)
{
tmp = x[i*n+k];
for (j = i + 1; j < n; j++)
tmp -= x[j*n+k] * a[i*n+j];
x[i*n+k] = tmp/a[i*n+i];
}
return 1;
}