LDLT-разложение

[nabbla1]

В некотором смысле, "заметки для себя": алгоритмы линейной алгебры обладают нехорошим свойством write-only, особенно если там участвует "треугольная адресация". Пока разбираешь - всё вроде понятно, а чуть-чуть погодя смотришь на это безобразие - и оно превращается в набор вложенных циклов и десяток индексов массива, которые обновляются на каждой

Comments

[quartisz]

А где можно найти готовую процедуру для LDLt на Pascal или Фортран, ну или на Си на худой конец?

[nabbla1]

У меня есть:

procedure CholeskyLDL (const a: array of Real; N: Integer; out u: array of Real);
var i: Integer; //счётчик для самого глубокого цикла
icol: Integer; //столбец, на котором текущий обраб. элемент
irow: Integer; //строка на которой текущий обраб. элемент
j: Integer; //указывает на элемент u[1, icol]
k: Integer; //указывает на элемент u[irow, icol], который и обрабатывается
l: Integer; //указывает на элемент u[i, irow]
m: Integer; //указывает на элемент u[i, icol]
d: Integer; //указывает на диагональный элемент u[i,i]
w: Real; //промежут. значение

begin
if N <= 0 then
raise Exception.Create('CholeskyLDL: N <= 0');

if ( Length(a) < N * ( N + 1 ) div 2 ) or (Length(u) < N * (N + 1) div 2) then

[quartisz]

Спасибо. А массивы почему одномерные?
Я вот такую написал процедуру.

type
TExtVec = array of Extended;
TExtMatr = array of TExtVec;

procedure LDLT(n: Integer; var a,L: TExtMatr; var d: TExtVec);
var i,j,k: Integer;
s: Extended;
begin
SetLength(L,n+1,n+1);
SetLength(d,n+1);
s:= 0;
for i:= 1 to n do
for j:= i+1 to n do
begin
s:= A[j][i];
for k:= 1 to i do
s:= s-L[i,k]*D[k]* L[j,k];
if (i = j) then
begin
D[i]:= s;//диагональный элемент
L[i,i]:= 1;//диагональ
end
else L[j,i]:= s/D[i];//внедиагональный элемент
end;
end;

[nabbla1]

Я брал пример с этого дела: https://people.sc.fsu.edu/~jburkardt/f77_src/asa007/asa007.html

когда всё это писалось, динамические массивы были редкостью, предпочитали всю эту конструкцию передать "одним куском", и в плане производительности это хорошо - заведомо всё рядом лежит, помещается сразу в кэш. Но адресация там ужасающая конечно...

[quartisz]

Сейчас и двумерные динамические массивы могут быть треугольными по памяти. Спасибо.

[quartisz]

Исправлена ошибка

procedure LDLT(n: Integer; var a,L: TExtMatr; var d: TExtVec);
var i,j,k: Integer;
s: Extended;
begin
SetLength(L,n+1,n+1);
SetLength(d,n+1);
s:= 0;
for i:= 1 to n do
for j:= i to n do
begin
s:= A[j][i];
for k:= 1 to i do
s:= s-L[i,k]*D[k]* L[j,k];
if (i = j) then
begin
D[i]:= s;
L[i,i]:= 1;
end
else L[j,i]:= s/D[i];
end;
end;