Сложность метода отражений для решения СЛАУ

[PhilB]

Добрый день.
Возникли некоторые сложности с реализацией или пониманием метода отражений для решений СЛАУ.
Реализую так.
1. Осуществляю QR-разложение.
На k-этапе этапе:
1) строю матрицу Хасухолдера.
H_k = E - 2w_kw_k^T,
где
w_k = \mu_k(a_{1k} - \beta_k; a_{2k}; ...; a_{nk})^T, \beta_k = sgn(-a_{1k}) \sqrt{\sum^n_{i=1} a_{ik}^2}, \mu_k = \frac{1}{2 \beta_k^2 - 2 \beta_k a_{kk}}
2) далее оперирую k-подматрицей (в которой вычеркнуты k строк и столбцов), получаю подматрицу A_k = H_kA, в которой обнулен 1 столбец подматрицы (без первого элемента).
2. Полученную на первом шаге треугольную матрицу использую для обратного хода.

Код функции:

Используется синтаксис C

double * solveSluUsingReflection(Slu slu) {

    size_t i,j,m,k;
    size_t n = slu.n;
    double ** a = malloc((n)*sizeof(double*));
    for (i = 0; i < n; i++) {
        a[i] = malloc((n+1)*sizeof(double));
        memcpy(a[i], slu.a[i], n*sizeof(double));
        a[i][n] = slu.b[i];
    }

    double ** h = getUnitMatrix(n);
    double ** e = getUnitMatrix(n);
    double * w = malloc(n * sizeof(double));
    double ** aa = allocMatrix(n,n);
    for (i = 0; i < n-1; i++) {
        double s = 0;
        for (k = i; k < n; k++) {
            s += a[k][i]*a[k][i];
        }
        if (fabs(s) < EPS) continue;
        double beta = (a[i][i] < 0) ? sqrt(s) : - sqrt(s);
        double mu = 1 / sqrt(2 * beta * beta - 2 * beta * a[i][i]);
        w[i] = mu * (a[i][i] - beta);
        for (j = i+1; j < n; j++) {
            w[j] = mu*a[j][i];
        }
        for (j = i; j < n; j++) {
            for (k = i; k < n; k++) {
                h[j][k] = e[j][k] - 2 * w[j] * w[k]; //ww[j][k];
            }
        }
        for (m = i; m < n; m++) {
            for (j = i; j < n+1; j++) {
                aa[m][j] = 0;
                for (k = i; k < n; k++) {
                    aa[m][j] += h[m][k] * a[k][j];
                }
            }
        }
        for (j = i; j < n; j++) {
            for (k = i; k < n+1; k++) {
                a[j][k] = aa[j][k];
            }
        }
    }

    // обратный ход
    double s;
    double * c = (double *)malloc(n * sizeof(double));
    for (i = n-1; i < n; i--) {
        s = a[i][n]; // b[i]
        for (j = i+1; j < n; j++) {
            s -= a[i][j] * c[j];
        }
        c[i] = s / a[i][i];
    }

    printf("\n");
    return c;

}
 

И вроде бы все отлично, но: получается сложность O(n^4); для СЛУ с 500 ур-ями, этот метод работает около 5 минут, в отличие от Гаусса, который работает около 1 секунды.
В книге Бахвалова и Амосова речь идет о O(n^3), и вообще, о том, что этот метод один из самых быстрых. В чем моя ошибка? Заранее благодарен.