
%  semidefinite_inverse_power.m   ver 1..0  by Tom Irvine

function[fn,Xhat]=semidefinite_inverse_power(M,K)

rng(1)
num=100;

%

    disp(' ');
    iv=1;
    K=K+iv*M;

%
n=max(size(M));
%

Y=zeros(n,1);
%
for i=1:n
    Y(i)=rand;
end
Y=Y-0.5;
Yb=Y;
%
Kinv=pinv(K);
%
for i=1:num
    Xhat=Kinv*Y;
    Yhat=M*Xhat;
    rho=Xhat'*Yhat/(Xhat'*Y);
    Y=Yhat/sqrt(Xhat'*Y);
    
    d=max(abs(Y-Yb));

    if(i>2 && d<1.0e-08)
        break;
    end    
    Yb=Y;
end
%
rho=1/rho;
%
rho=rho-iv;
%

if(abs(rho)< 1.0e-04)
    rho=0.;
end

%
out1=sprintf(' The minimum eigenvalue = %8.4g ',rho);
disp(out1);
out1=sprintf('              frequency = %8.4g Hz',sqrt(abs(rho))/(2*pi));
disp(out1);
disp(' ');

    r=Xhat'*M*Xhat;
    Xhat=Xhat/sqrt(r);


if(Xhat(1)<0)
    Xhat=-Xhat;
end  

fn=sqrt(rho)/(2*pi);
