
% Example of memory corruption in Octave's ichol function.
% This script constructs a symmetric positive definite matrix A
% as a sparse matrix and proceeds to compute its Cholesky factor L.
% The Cholesky factorisation is done using ichol with type "ict"
% (incomplete Cholesky with thresholding). The ichol builtin dispatches to a
% C++ fuction named "icholt".
% A zero drop tolerance essentially means that the incomplete factorization
% conicides with the exact Cholesky factorization.
% Octave 4.2.1 crashes when running this script in Linux, resulting in one of
% the following:
%   - Command window freezes, does not accept input,
%   - an error message saying "Memory corruption" occurs on the console,
%   - Octave aborts with "SIGNAL 6", i.e. SIGABRT under Linux.
% See ichol_mark_array.patch for a modification of the icholt function in file
%   __ichol__.cc
% from the Octave source tree for a possible way how this problem can be
% fixed. After applying this patch to the mentioned source file this example
% runs without any noticable problems giving a very small norm(A-L*L.') as
% expected.
% The problem occurs when icholt computes the third column of the Cholesky factor.
% The third row of L is full when the processing of the third row begins.
% After substracting the first row from the third, all elements in the third
% row are eliminated due to the construction of A.
% When the second row is constructed from the third the third row gets
% non-zero again. The corresponding column indices are inserted to the "vec"
% array a second time. This results in about 2 * n indices in "vec" which is
% of length only n.

k=4;
n=2^4;
Afull=diag(ones(n,1))/2+triu([zeros(n,2),[n.^-[1;1;2]*ones(1,n-2);zeros(n-3,n-2)]],1);
A=sparse(Afull+Afull.');
opts.type="ict";
opts.droptol=0;
L=ichol(A,opts);
norm(A-L*L.')
