function [lamda,kpc]=pdsolv(nli,mwi,sf)
%PDSOLV is a utility function for FPSEUDO
%       [lamda,Kpc]=PDSOLV(NLI,MWI,sf) solves the real symmetric eigenvalue
%       problem  (MWI.*sf)*Kpc = lamda*(NLI.*sf)*Kpc  using singular value
%       decompositions.  When NLI is almost singular small changes in the
%       scalar sf affect the result.
%       Kpc a vector.

%       Dr M.P. Ford 2nd September 1987
% Copyright (c) 1987 by GEC Engineering Research Centre & Cambridge Control Ltd
%       MRN0019

[x,d,temp]=svd(nli.*sf);   % since N is real symmetric
			   % singular values = eigenvalues
			   % but SVD has better numerical properties.
if any(diag(d)==0);   % a perfect fit
   crlf=setstr([13,10]);
   error(['Matrix is singular to working precision', crlf,...
	  'Order of compensator too high or not enough frequency points?'])
end
id=diag(sqrt(diag(d)).\1);   % inverse of square roots of d  (diagonal)
x=real(x);             % x should be real
mwid=id*x'*(mwi.*sf)*x*id;
mwid=real((mwid+mwid')./2);  % make sure it is real symmetric
[s,lamda,temp]=svd(mwid);
[lamda,i]=min(diag(lamda));  % last one should be the minimum
kpc=x*id*real(s(:,i)); %  s should be real also