// You are free to use, redistribute and modifiy this code if you include
// as a comment the author of the original code 
// (c) Sandor Zoltan Nemeth, 2009.
function p=piso(x,E)
   [n,n]=size(E)
   if det(E)==0 then 
      error ('The input cone must be an isotone projection cone'); return
   else
   V=inv(E) 
   U=-V'
   F=-V*U
   G=F-diag(diag(F))
      for i=1:n
         for j=1:n
            if G(i,j)>0 then
               error('The input cone must be an isotone projection cone'); return
            end
         end
      end
      I=[1:n]
      n1=n-1
      cont=1
      for k=0:n1
        //disp(cont,"cont=")
        [q1,l]=size(I)
	     E1=E
	     I1=I
	     if l-1>=1
	       highest=I(l)
	       if highest<n
	          for h=n:-1:highest+1
		    E1(:,h)=zeros(n,1)
	          end 
	       end
               for j=l-1:-1:1 
                  low=I(j)+1
                  high=I(j+1)-1
                  if high>=low
                     for m=high:-1:low
                     E1(:,m)=zeros(n,1)

                  end
               end
	       lowest=I(1)
	       if lowest>1
	          for w=lowest-1:-1:1
	             E1(:,w)=zeros(n,1)
 	          end
	       end
             end
       end
       if l==1 
	E1=zeros(n,n)
	E1(:,I(1))=E(:,I(1))
       end
     
         V1=pinv(E1)
         U1=-V1'
       //disp(x,"start x=")
       //disp(I,"I=")
   
	  //disp(U1,"U1=")
         //disp(V1,"V1=")
         for j=l:-1:1
          zz=x'*U1(:,I(j))
	  if zz>0 I1(j)=[]
            end       
         end
       //disp(I1,"I1=")
         [q2,ll]=size(I1)
            if cont==0 then p=x; return
            elseif ll==0 then p=zeros(n,1); return 
            else 
         //disp(E1,"E1=")
         //disp(V1,"V1=")
               A=E1*V1
         //disp(A,"A=")
               x=A*x
               //disp(x,"x=")
               p=x;
            end
            if ll==l cont=0
            else cont=1
            end 
         I=I1
      end
   end
endfunction
