function [h,InfResp] = fresp(D,w)
% Frequency response of SPARSS model.
% InfResp is a logical scalar flagging singularities in the frequency
% response.

%	 Author(s): P.Gahinet 
%   Copyright 2020 The MathWorks, Inc.

% Form vector s of complex frequencies
[s,sn] = ltipack.utGetComplexFrequencies(w,D.Ts);

% Seek arrow structure
D = xsort(D);

% Scaling
if D.Ts==0
   % Scale based on mid-range frequency
   f = sort(abs(w(w>0 & isfinite(w))));
   if ~isempty(f)
      fmid = sqrt(f(1)*f(end));
      D = scale(D,[],abs(D.a)+fmid*abs(getE(D)));
   end
else
   D = scale(D,[],abs(D.a)+abs(getE(D)));
end

% Extract data
a = D.a;
b = D.b;
c = D.c;
d = full(D.d);
e = getE(D);
nf = numel(w);
[ny,nu] = size(d);

% Configure solver and parallelism
OPTS = getSolverConfig(D);
UseParallel = ltipack.util.checkParallel(OPTS.UseParallel);
if UseParallel
   ParPool = gcp;
else
   ParPool = matlab.internal.threadPool;
end

% Performance optimization to minimize number of backsolves
if ny<nu
   cc = full(c)';
else
   cc = c';  b = full(b);
end

% Compute response
h = zeros(ny,nu,nf,'like',1i);
parfor (ct=1:nf, ParPool)
   h(:,:,ct) = localFRESP(a,b,cc,d,e,s(ct));
end
if allfinite(s)
   InfResp = ~allfinite(h);
else
   InfResp = ~allfinite(h(:,:,isfinite(s)));
end

% Factor in delays
Delay = D.Delay;
nfd = numel(Delay.Internal);
if nfd>0
   ny = ny-nfd;  nu = nu-nfd;
   h = ltipack.util.addInternalDelayResp(h,Delay.Internal,sn);
end
if any(Delay.Input) || any(Delay.Output)
   Tdio = Delay.Input(:,ones(1,ny)).' + Delay.Output(:,ones(1,nu));
   h = ltipack.util.addIODelayResp(h,Tdio,sn);
end      
   
%--------------------------------------------
function h = localFRESP(a,b,cc,d,e,s)
% Computes frequency response D+C*(sE-A)\B. Here CC=C'.
[ny,nu] = size(d);
if isfinite(s)
   % Note:
   % 1) Evaluating c*x is most efficient when c is sparse and x full
   % 2) c'*x is faster than c*x for sparse c
   if ny<nu
      % Evaluate as D+(B'*((s*E-A)'\CC))'
      x = matlab.internal.math.nowarn.mldivide(s*e-a,cc,true);
      if allfinite(x)
         h = d + (b'*x)';
      else
         h = inf(ny,nu);
      end
   else
      % Evaluate as D+CC'*((s*E-A)\B)
      x = matlab.internal.math.nowarn.mldivide(s*e-a,b,false);
      if allfinite(x)
         h = d + cc'*x;
      else
         h = inf(ny,nu);
      end
   end      
else
   h = nan(ny,nu);
end
