LDR разложение, алгоритм.

[shukshin]

Добрый день, уважаемые участники форума!
Не могли бы вы мне подсказать (или дать источник, где это есть) про алгоритм разложения матрицы в произведение LDR (нижнетреугольная, диагональная, верхнетреугольная) для решений СЛАУ. Сам в интернете нигде не найду.
Материалы из линейной алгебры дают очень запутанный и бессмысленный рекурсивный алгоритм. (Хотя бы тем, что для вычислений там необходимо вычислять обратную матрицу.)

[мат-ламер]

shukshin в сообщении #907565 писал(а):

Сам в интернете нигде не найду.

Учебники не пробовали смотреть (Воеводин?)?

[shukshin]

там есть на эту тему, но все-таки не оно (разложение $A = L L^T $ вроде)

[arseniiv]

Обычно это, насколько я понимаю, называется LDU — может, потому плохо искалось?

[Yuri Gendelman]

Уилкинсон, Райнш. Справочник алгоритмов на языке АЛГОЛ. Линейная алгебра.
Алгоритм I.7 на стр.92 (разложение LDU).

[shukshin]

Нашел вот здесь в итоге:


Но не получается алгоритм, ни в какую. Может свежим взглядом, кто зацепит? (индексы сдвинуты на 1, так как в С нумерация массивов с 0)

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

void LDU_decomposition(matrix_t *A, matrix_t *L, matrix_t *D, matrix_t *U)
{
        int i, j ,k;
        double S;
 
        Make_E(L), Make_E(D), Make_E(U);
 
        D->p[0][0] = A->p[0][0];
 
        for(j = 1; j < A->x; ++j)
        {
        
          for(i = 0; i <= j - 1; ++i)
          {
                    if(i == 0)
                        {
                                U->p[i][j] = (double)(A->p[i][j]) / (D->p[i][i]);
                                L->p[j][i] = (double)(A->p[j][i]) / (D->p[i][i]);
                        }
                        else if(i > 0)
                        { 
                                for(S = 0, k = 0; k < i - 1; ++k)
                                        S += L->p[i][k] * D->p[k][k] * U->p[k][j];
                                U->p[i][j] = (double)(A->p[i][j] - S) / (D->p[i][i]);
 
                                for(S = 0, k = 0; k < j - 1; ++k)
                                        S += U->p[k][i] * D->p[k][k] * L->p[j][k];
                        L->p[j][i] = (double)(A->p[j][i] - S) / (D->p[i][i]);
                        }
          }
         
          for(S = 0, k = 0; k < j - 1; ++k)
                S += (U->p[k][j]) * (D->p[k][k]) * (L->p[j][k]);
          D->p[j][j] = (A->p[j][j] - S);
          
         
        }
}

[мат-ламер]

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