function [F,di] = factorise(C,SVDONLY)
% factorise  Use Cholesky or SVD to factorise covariance matrices.
%
% Backend IRIS function.
% No help provided.

% -IRIS Toolbox.
% -Copyright (c) 2007-2012 Jaromir Benes.

try
    SVDONLY; %#ok<VUNUS>
catch %#ok<CTCH>
    SVDONLY = false;
end

%**************************************************************************

csize = size(C);
C = C(:,:,:);
nx = size(C,1);
nloop = size(C,3);
F = nan(size(C));
di = nan(1,nloop);

for i = 1 : nloop
   if i > 1 && all(all(C(:,:,i) == C(:,:,i-1)))
      % Copy the previous factor matrix if the cov matrix is identical to
      % the previous cov matrix.
      F(:,:,i) = F(:,:,i-1);
      if nargout > 1
         di(i) = di(i-1);
      end
      continue
   end
   Ci = (C(:,:,i)+C(:,:,i).')/2;
   cholFailed = false;
   if ~SVDONLY
      % Unless declined by the user, attempt Cholesky first. Cholesky is
      % faster than SVD but may fail because of the matrix being
      % numerically non-positive-definite.
      try
         F(:,:,i) = chol(Ci).';
      catch %#ok<CTCH>
         cholFailed = true;
      end
   end
   if SVDONLY || cholFailed
      % If requested by the user or if Cholesky failes, use SVD which works
      % under any circumstances. Set negative svd values
      [U,S] = svd(Ci);
      S = diag(S).';
      S(S < 0) = 0;
      % Impose sign convention that first element of each singular vector
      % is positive. This is to make results identical across different
      % processors.
      index = sign(U(1,:)) == -1;
      U(:,index) = -U(:,index);
      S = sqrt(S);
      F(:,:,i) = U .* S(ones([1,nx]),:);
   end
   if nargout > 1
      di(i) = maxabs(Ci - F(:,:,i)*F(:,:,i).');
   end
end

if length(csize) > 3
   F = reshape(F,csize);
end

end
