Разложение методом Холецкого

Методом Холецкого.
function [A,p] = chol(A)
% Embedded MATLAB Library function.
%
% Limitations:
% 1) Does not support CHOL with two output arguments
% 2) May factor matrices that MATLAB will not.
% 3) May fail to factor matrices that MATLAB will factor.

% $INCLUDE(DOC) toolbox\matlab\matfun\chol.m$
% Copyright 2002-2004 The MathWorks, Inc.
% $Revision: 1.1.4.2 $ $Date: 2004/07/01 04:50:25 $

% Translated from LAPACK functions zpotf2.f/cpotf2.f/dpotf2.f/spotf2.f

eml_assert(isfloat(A),['Function ''chol'' is not defined for values of class ''' class(A) '''.']);
eml_assert(nargin == 1, 'Not enough input arguments.');
eml_assert(size(A,1) == size(A,2), 'Matrix must be square.');
eml_assert(nargout<2,'Embedded MATLAB''s CHOL function does not support the two output version.');

if isempty(A)
return;
end

n = size(A,1);
A = triu(A);
p = 0;
for j = 1 : n

%s = real(A(j,j)) - A(1:j-1,j)' * A(1:j-1,j);
t = cast(0,class(A));
if j>1
for k = 1 : j-1
akj = A(k,j);
% t = t + real(conj(A(k,j)) * A(k,j));
t = t + real(akj)*real(akj) + imag(akj)*imag(akj);
end
end
ajj = real(A(j,j)) - t;

if ajj <= 0
A(j,j) = ajj;
p=j;
if nargout <= 1
error(['Not positive definite. For a matrix with eigenvalues ' ...
'close to EPS, MATLAB and Embedded MATLAB may ' ...
'disagree on whether or not it is positive definite ' ...
'due to rounding errors.']);
end
return;
end

ajj = sqrt(ajj);
A(j,j) = ajj;

if j < n
%A(j,j+1:n) = A(j,j+1:n) - A(1:j-1,j)'*A(1:j-1,j+1:n);
for k = j+1:n
if isreal(A)
temp = cast(0,class(A));
if j > 1
for w = 1 : j-1
temp = temp + A(w,j) * A(w,k);
end
end
A(j,k) = (A(j,k) - temp)/ajj;
else
tre = cast(0,class(A));
tim = cast(0,class(A));
tim2 = cast(0,class(A));
if j > 1
for w = 1 : j-1
awj = A(w,j);
awk = A(w,k);
tre = tre + real(awj) * real(awk) + imag(awj) * imag(awk);
tim = tim + real(awj) * imag(awk);
tim2 = tim2 + imag(awj) * real(awk);
end
end
A(j,k) = (A(j,k) - complex(tre,tim - tim2))/ajj;
end
end
end
end

function B = cast(A,cls)

switch cls
case 'double',
B = double(A);
case 'single',
B = single(A);
otherwise,
eml_assert(false,'error','Internal error: unexpected class.');
end
Только это не Паскаль!!!!!! Это стандартная функция из MATLAB.