Не работает алгоритм QR разложения матриц c++

Вот мой код:

#include<iostream>
#include<fstream>
#include<vector>
#include<cmath>


using namespace std;


vector < vector < double>> multiply(vector < vector < double>> first, vector < vector < double>> second, int N)
{
  vector < vector < double>> res(N, vector<double>(N));
  for (int i = 0; i < N; i++)
  {
    for (int j = 0; j < N; j++)
    {
      for (int k = 0; k < N; k++)
      {
        res[i][j] += first[i][k] * second[k][j];
      }
    }
  }
  return res;
}


vector < vector < double>> multiply(vector<double> first, vector<double> second, int N)
{
  vector < vector < double>> res(N, vector<double>(N));
  for (int i = 0; i < N; i++)
  {
    for (int j = 0; j < N; j++)
    {
      res[i][j] = first[i] * second[j];
    }
  }
  return res;
}


int sign(double a)
{
  if (a > 0)
  {
    return 1;
  }
  else if (a < 0)
  {
    return -1;
  }
  else
  {
    return 0;
  }
}


int main()
{
  int dimension;
  {
    ifstream fin("dimension.txt");
    fin >> dimension;
  }
  vector < vector < double>> matrix(dimension, vector<double>(dimension)), Q(dimension, vector<double>(dimension, 0)), R(dimension, vector<double>(dimension)), H(dimension, vector<double>(dimension)), E(dimension, vector<double>(dimension, 0));
  {
    ifstream fin("matrix.txt");
    for (int i = 0; i < dimension; i++)
    {
      for (int j = 0; j < dimension; j++)
      {
        fin >> matrix[i][j];
      }
    }
  }
  int accuracy;
  {
    ifstream fin("accuracy.txt");
    fin >> accuracy;
  }
  vector<double> V(dimension);


  for (int i = 0; i < dimension; i++)
  {
    E[i][i] = 1;
    Q[i][i] = 1;
  }


  //алгоритм QR
  double sum = -1;
  while (accuracy >= sum)
  {
    for (int i = 0; i < dimension - 1; i++)
    {
      //определение компонент вектора V
      for (int j = 0; j < dimension; j++)
      {
        if (j < i)
        {
          V[j] = 0;
        }
        else if (j == i)
        {
          double evkl = 0;
          for (int k = i; k < dimension; k++)
          {
            evkl += pow(matrix[k][i], 2);
          }
          evkl = sqrt(evkl);
          V[j] = matrix[j][i] + sign(matrix[j][i]) * evkl;
        }
        else
        {
          V[j] = matrix[j][i];
        }
      }

      //получаем Н
      double VtV = 0;
      for (int i = 0; i < dimension; i++)
      {
        VtV += V[i] * V[i];
      }
      H = multiply(V, V, dimension);
      for (int i = 0; i < dimension; i++)
      {
        for (int j = 0; j < dimension; j++)
        {
          H[i][j] = E[i][j] - 2 * H[i][j] / VtV;
        }
      }

      //получаем A(matrix) и Q 
      //начиная с этого момента вычисления становятся некорректными
      matrix = multiply(H, matrix, dimension);
      Q = multiply(Q, H, dimension);
    }

    sum = 0;
    matrix = multiply(matrix, Q, dimension);
    for (int i = 0; i < dimension - 1; i++)
    {
      for (int j = i + 1; j < dimension; j++)
      {
        sum += pow(matrix[j][i], 2);
      }
    }
    sum = sqrt(sum);

  }


}

По какой-то причине неверно перемножаются матрицы matrix и H в функции multiply(после перемножения поддиагональные элементы матрицы должны обнулятся, но этого не происходит), хотя функция перемножения написана правильно, я её уже неоднократно использовал. Все вычисления до этого этапа верные, их я проверял. Помогите понять, в чём заключается ошибка.