function R = gram(D,type,Options)
% Computes Cholesky factor of
%    * Controllability Gramian if TYPE = 'c' or 'cf'
%    * Observability Gramian if TYPE = 'o' or 'of'
%
%   Note: Third input OPTIONS must be of class ltioptions.gram.

%   Author(s): P. Gahinet, S. Gumussoy
%   Copyright 1986-2021 The MathWorks, Inc.
Ts = D.Ts;
if hasInternalDelay(D)
   throw(ltipack.utNoDelaySupport('gram',Ts,'internal'))
end

if isempty(Options.TimeIntervals) && isempty(Options.FreqIntervals)
   % Standard gramian computation

   % Check that E is nonsingular
   [isProper,Dr] = isproper(D,false);
   if ~isProper || order(Dr)<order(D)
      error(message('Control:general:NotSupportedSingularE','gram'))
   end
   
   % Descriptor only supported for real data
   if isreal(D)
      % Use available form
      [a,b,c,~,~,e] = getABCDE(D);
   else
      % Derive explicit form
      [a,b,c] = getABCD(D);   e = [];
   end
   
   % Solve Lyapunov equation
   if Ts==0
      LyapSolver = @lyapchol;
   else
      LyapSolver = @dlyapchol;
   end
   try
      if any(type=='c')
         R = LyapSolver(a,b,e);
      else
         R = LyapSolver(a',c',e');
         if ~isempty(e)
            % Wo = (R*E)'*(R*E)
            R = R*e;
            if any(type=='f')
               [~,R] = qr(R);
            end
         end
      end
   catch %#ok<CTCH>
      %Some eigenvalue of A or A' lies in the unstable region.
      error(message('Control:transformation:gramerror'))
   end
else
   % Limited gramian computation

   % Make model explicit
   [isProper,Dr] = isproper(D,true);
   if ~isProper || order(Dr)<order(D)
      error(message('Control:general:NotSupportedSingularE','gram'))
   end
   
   if any(type=='c')
      R = gramLimited(Dr.a,Dr.b,Dr.Ts,Options);
   else
      R = gramLimited(Dr.a',Dr.c',Dr.Ts,Options);
   end
end

