function [v,l] = pdeeig(varargin) % b,p,e,t, varargin
%PDEEIG Solve PDE eigenvalue problem.
%   [V,L]=PDEEIG(MODEL,C,A,D,R) solves the FEM formulation of the PDE
%   eigenvalue problem -div(C*grad(U))+ A*U= L*D*U, using the geometry,
%   mesh, and boundary conditions described by the MODEL, a PDEMODEL
%   object.
% 
%   The coefficients C, A, and D of the PDE problem can be given in a wide
%   variety of ways. See 'PDE Coefficients' page of the documentation for
%   details.
% 
%   R is a two element vector, indicating an interval on the real axis.
%   (The left-hand side may be -Inf.) The algorithm returns all eigenvalues
%   in this interval in L.
% 
%   V is a matrix of eigenvectors. For a single equation eigenvalue
%   problem, each column in V is an eigenvector evaluated at NP nodes given
%   by MODEL.MESH.NODES. For a system of N equations with NP nodes, the
%   first NP rows of V corresponds to eigenvectors of first equation
%   evaluated at NP nodes, following NP rows corresponds to second equation,
%   and so on for N equations. Thus, for a system of equations there will be
%   NP*N number of rows in V.
% 
%   [V,L]=PDEEIG(B,P,E,T,C,A,D,R) allows the geometry and mesh to be
%   described by P, E, and T, see INITMESH for details. B describes the
%   boundary conditions of the PDE problem. Refer to the documentation page
%   'Boundary Conditions by Writing Functions' for details.
% 
%   [V,L]=PDEEIG(Kc,B,M,R) produces the solution to the generalized sparse
%   matrix eigenvalue problem Kc*UI=L*B'*M*B*UI, U=B*UI.
%
%   Note: This function does not support quadratic triangular elements
%

%       Copyright 1994-2017 The MathWorks, Inc.

if nargin > 0
    [varargin{:}] = convertStringsToChars(varargin{:});
end

narginchk(4,10);
% b,p,e,t, varargin
% [V,L]=PDEEIG(K,B,M,R)

b = varargin{1};
meshedpde=false;      
if isa(b, 'pde.PDEModel') 
    if isa(b.Mesh, 'pde.FEMesh')
      meshedpde=true;
      mypde = varargin{1};    
      msh = mypde.Mesh;
      [p,e,t] = msh.meshToPet(); 
      b = pde.PDEModel(mypde.PDESystemSize);  
      b.BoundaryConditions = mypde.BoundaryConditions;
      if nargin > 1
        varargin = {varargin{2:end}};
      else
        varargin = {};
      end
      quadraticTriNotSupported(t);
    elseif ~b.IsTwoD
      error(message('pde:parabolic:pdemodelNoMesh'));  
    end   
end

parser = inputParser;
isCoeffForm = ~issparse(b);
if(isCoeffForm)
 if ~meshedpde
  parser.addRequired('b', @pdeIsValid.b);
  parser.addRequired('p', @pdeIsValid.p2d);
  p = varargin{2};
  e = varargin{3};
  eTest = @(e) pdeIsValid.meshGeomAssoc(e, p);
  parser.addRequired('e', eTest);
  parser.addRequired('t', @pdeIsValid.t2d);
 end
  parser.addRequired('c', @pdeIsValid.coefficient);
  parser.addRequired('a', @pdeIsValid.coefficient);
  parser.addRequired('d', @pdeIsValid.coefficient);
  parser.addRequired('r', @isnumeric);
  parser.addParameter('Stats', 'on', @pdeIsValid.onOrOff);  
  parser.parse(varargin{:});
  pres = parser.Results; 
  c = pres.c;
  a = pres.a;
  d = pres.d;
  r = pres.r;
  if ~meshedpde  
      b = pres.b;
      p = pres.p;
      e = pres.e;
      t = pres.t;
      quadraticTriNotSupported(t);
  end    
  np=size(p,2);
  % Boundary contributions
  [Q,G,H,R]=assemb(b,p,e);
  % Number of variables
  N=size(Q,2)/np;
  if(isa(b, 'pde.PDEModel'))
    thePde = b;
  else
    thePde = pde.PDEModel(N);
  end
  [K,M,F] = thePde.assema(p,t,c,a,zeros(N,1));
  [K,~,B,~] = assempde(K,M,F,Q,G,H,R);
  [~,M] = thePde.assema(p,t,0,d,zeros(N,1));
else
  parser.addRequired('K', @issparse);
  parser.addRequired('B', @issparse);
  parser.addRequired('M', @issparse);
  parser.addRequired('r', @isnumeric);
  parser.addParameter('Stats', 'on', @pdeIsValid.onOrOff);
  parser.parse(varargin{:});
  pres = parser.Results; 
  K = pres.K;
  B = pres.B;
  M = pres.M;
  r = pres.r;  
end
stats = pres.Stats;

M=B'*M*B;

if (issymmetric(K)&&issymmetric(M)&&r(1)~=-Inf) 
    % Positive definite is checked inside sptarn
    spd=1;
else
    spd=0;
end

[v,l,ires]=sptarn(K,M,r(1),r(2),spd, 'Stats', stats);

if ires<0
  error(message('pde:pdeeig:MoreEigenvaluesMayExist'));
end

if ~isempty(v)
  v=B*v;
end

end

