%MKLDL  Generates code for LDL inversion
% F = MKLDL(N) Prototype
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
function f = mkldl(n)
 if nargin == 0
  n = 6;
  a = magic(n);
  a = a'*a; % a = a(a>0);
  [d,id] = deal(zeros(1,n));
  [l,w]  = deal(eye(n));
  f = mkldl(n)
  eval(f)
  norm(l*diag(1./id)*l'-a)
  return
 end
 mode = 1;
 [i,j,k] = deal(0);
 f = newline;
 for j = 1:n
  f = rep(['% J=j' newline]);
  if j == 1
   f = rep(['id(j) = 1/a(j,j);' newline]);
  else
   f = rep('id(j) = a(j,j)');
   for k = 1:j-1
    if k == 1
     f = rep(' - a(j,k)*l(j,k)');
    else
     f = rep(' - w(jw,kw)*l(j,k)');
    end
   end
   f = rep([';' newline 'id(j) = 1/id(j);' newline]);
  end
  for i = j+1:n
   f = rep(['% I=i' newline]);
   if j == 1
    f = rep(['l(i,j) = a(i,j)*id(j);' newline]);
   else
    f = rep('w(iw,jw) = a(i,j)');
    for k = 1:j-1
     if k == 1
      f = rep(' - l(i,k)*a(j,k)');
     else
      f = rep(' - l(j,k)*w(iw,kw)');
     end
    end
    f = rep([';' newline 'l(i,j) = w(iw,jw)*id(j);' newline]);
   end
  end
 end
 
 function s = rep(s)
  if mode == 1
   s = regexprep(s,...
    {'\<i\>,\<j\>' '\<j\>,\<j\>' '\<i\>,\<k\>' '\<j\>,\<k\>'},...
    {s2(i,j,n) s2(j,j,n) s2(i,k,n) s2(j,k,n)});
   s = regexprep(s,...
    {'\<iw\>,\<jw\>' '\<jw\>,\<jw\>' '\<iw\>,\<kw\>' '\<jw\>,\<kw\>'},...
    {s2(i-2,j-1,n-2) s2(j-2,j-1,n-2) s2(i-2,k-1,n-2) s2(j-2,k-1,n-2)});
  end
  s = [f regexprep(s,{'\<iw?\>' '\<jw?\>' '\<kw?\>'},{int2str(i) int2str(j) int2str(k)})];
 end
 
 function s = s2(a,b,n)
  if a < b, 
   s = s2(b,a,n);
  else
  s = int2str(a+(2*n-b)*(b-1)/2);
  end
 end
end