Прямой и обратный ход Гаусса в одном цикле

#cpp #алгоритм

Вопрос больше теоретического характера. Каким образом можно объединить прямой и обратный
ход Гаусса при решении СЛАУ в одну процедуру, то есть имеется в виду, в один цикл,
в один "проход"?

Еще раз повторюсь, мне не нужен целиком алгоритм. Мне нужна лишь идея на словах с
очень кратким пояснением к ней.

Быть может, попробовать организовать какой-нибудь цикл от n-1 до 0 (ну если СЛАУ
nxn), где магическим образом будут и коэффициенты матрицы зануляться, и вектор решения
считаться? 

Если вдруг кому может понадобиться кусочек кода, реализующий эти два хода, то оно вот:

/// Прямой ход
bool gausss( tdouble ** mtx, tdouble * vect )
{
    tdouble c;
    for ( int j = 0; j < n; j++ ) {
        tdouble temp = abs( mtx[j][j] );
        int mem = j; максимум
        for ( int i = j + 1; i < n; i++ ) { // от j+1
            if ( temp < abs( mtx[i][j] ) ) {
                temp = abs( mtx[i][j] );
                mem = i;
            }
        }

        if ( temp < eps ) return false;

        changing( mtx, vect, mem, j );

        for ( int i = j + 1; i < n; i++ ) {
            c = mtx[i][j] / mtx[j][j];
            for ( int k = j; k < n; k++ ) {
                mtx[i][k] -= c * mtx[j][k];
            }
            vect[i] -= c * vect[j];
        }
    }
    return true;
}

/// Обратный ход
void gausss_reverse( tdouble ** mtx, tdouble * vect, tdouble * sol )
{
    for ( int i = n - 1; i >= 0; i-- ) {
        tdouble temp = 0;
        for ( int j = i + 1; j < n; j++ ) {
            temp += mtx[i][j] * sol[j];
        }
        sol[i] = (vect[i] - temp) / mtx[i][i];
    }
}

    

Ответы

Ответ 1

Ну смотрите.

В чём заключается обычно метод Гаусса? Вы на прямом ходе получаете нули под диагональю.
А на обратном — над диагональю. Это можно легко объединить.

В объединённом алгоритме вы сначала выполняете шаг прямого метода — при этом под
диагональю получаются нули. Теперь вы можете текущей строкой «убить» ненулевые коэффициенты
над ней. Таким образом, у вас матрица будет выглядеть так (пример для третьего шага):

a11  0   0  a14
 0  a22  0  a24
 0   0  a33 a34
 0   0   0  a44
 0   0   0  a54
 ..............


На каждом шаге у вас добавляется один столбец в диагональный вид матрицы.

В принципе, вы можете на каждом шаге ещё и делить на диагональный элемент, чтобы
получить там единицу. Тогда в конце у вас столбец свободных членов станет столбцом-решением.