/usr/share/euler/progs/choleski.e is in euler 1.61.0-10.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 | comment
Choleski decomposition and LGS solver.
Eigenvalue iteration using Choleski decomposition
choleski(A)
lsolve(L,b)
choleigen(A)
endcomment
function choleski (A)
## Decompose a real matrix A, such that A=L.L'. Returns L.
## A must be a positive definite symmetric and at least 2x2 matrix.
n=cols(A);
L=zeros(size(A));
L[1,1]=sqrt(A[1,1]); L[2,1]=A[1,2]/L[1,1];
L[2,2]=sqrt(A[2,2]-L[2,1]^2);
loop 3 to n;
c=flipy(lusolve(flipx(flipy(L[1:#-1,1:#-1])),flipy(A[1:#-1,#])));
L[#,1:#-1]=c';
L[#,#]=sqrt(A[#,#]-c'.c);
end;
return L;
endfunction
function lsolve (L,b)
## Solve L.L'=b
return lusolve(L',flipy(lusolve(flipx(flipy(L)),flipy(b))));
endfunction
function choleigen (A)
## Iterates the Choleski-iteration, until the diagonal converges.
## A must be positive definite symmetric and at least 2x2.
if isvar("eps"); localepsilon(eps); endif;
L=choleski(A);
B=L'.L;
d=diag(B,0);
repeat
L=choleski(B);
B=L'.L;
dnew=diag(B,0);
if dnew~=d; break; endif;
d=dnew;
end;
return dnew;
endfunction
|