function [r]=ptest(a)

% Return r=1 if `a' is a real P-matrix (r=0 otherwise).  

if ~isreal(a) 
   'Matrix is not real', break
end 
n=length(a);
if a(1,1)<=0
     r=0;
elseif n==1
     r=1;   
else
     for i=2:n
           d(i,1)=a(i,1)/a(1,1);
           for j=2:n
              b(i-1,j-1)=a(i,j);
              c(i-1,j-1)=a(i,j)-d(i,1)*a(1,j);
           end
     end  
     r1=ptest(b);
     r2=ptest(c);
     r=r1*r2;
end