В чем заключается этот метод?
Умножение матрицы

на некоторый вектор

сводится к тому, что из этого вектора вычитается его удвоенная проекция на направление вектора

или, что то же самое, вектор

зеркально отражается относительно плоскости, перпендикулярной вектору

. При этом получается некоторый вектор

, длина которого равна длине

. Для каждой пары векторов

и

одинаковой длины вектор

определяется однозначно (с точностью до знака, но это не важно):

.
Системы уравнений здесь вот при чём. Основная трудность при их решении -- это приведение матрицы

системы к треугольному виду, когда ниже диагонали стоят нули (после этого элементы выше диагонали уничтожаются уже очень быстро).Так вот, на первом шаге подбирается такой вектор отражения

, который переводит первый столбец

матрицы

в столбец

, т.е. такой, у которого ниже диагонали стоят нули. Естественно, при этом

, а знак выбирается противоположным знаку первой компоненты столбца

(так, чтобы после вычитания

получить результат побольше; это уменьшает погрешности округления). На втором шаге аналогично подбираем отражение, уничтожающее элементы второго столбца матрицы (той, которая получилась после первого отражения), лежащие ниже диагонали; при этом в отражении участвуют только компоненты столбцов, начиная со второй, а первые не меняются. И т.д.