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

[nabbla1]

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

Comments

[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
raise Exception.Create('CholeskyLDL: not enough space in arrays');

j := 1;
k := 0;
// Factorize column by column, ICOL = column number.
for icol := 1 to N do begin
l := 0;

// IROW = row number within column ICOL.
for irow := 1 to icol do begin
inc(k);
w := a[k-1];
m := j;
d := 0;
for i := 1 to irow - 1 do begin
inc(l);
w := w - u[l-1] * u[m-1] * u[d];
inc(m);
inc(d, i+1);
end;
inc(l);

if irow = icol then
break;
if u[l - 1] <> 0.0 then
u[k - 1] := w / u[l - 1]
else begin
raise Exception.Create('CholeskyFactorization: matrix is not positive defined');
end;

end;

// End of row, estimate relative accuracy of diagonal element.
u[k - 1] := w;
inc(j, icol);
end;
end;

Здесь имеется в виду "треугольная адресация", т.е к примеру для матрицы 6х6 массивы a, u должны иметь длину 6*7/2 = 21, и там плотно записаны слева направо сверху вниз.

входная и выходные могут быть различными, а могут совпадать - вычисления будут коррентно выполнены "на месте".

[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;