Can you please help to find out why I am obtaining a negative PSNR value in the following Matlab code ?

Dear Experts,

I am getting a negative PSNR value with the following Matlab code. Can you please help ?

% main program to demonstrate how NEDI works
% Author: Xin Li
% Reference: New edge-directed interpolation 
% Li, X and Orchard, M T 
% IEEE Transactions on Image Processing. Vol. 10, no. 10, pp. 1521-1527. Oct. 2001 

clear all
close all

B = double(imread('RGB_IMAGE_X_2.JPG'));      % B is a color image of the size(B) = 256 x 256 
x = double(imread('RGB_IMAGE_X_1.JPG'));      % x is a color image of the size(x) = 128 x 128 

imshow(x/255,[]);
title('original LR image');

for i=1:3
    y(:,:,i)=sri(x(:,:,i),1);
end

figure, imshow(y/255,[]);
title('interpolated HR image');

[PSNR,SNR] = psnr(y,B)

PSNR =

-19.5022

SNR =

   17.9484

My question is why do I get the negative PSNR value ? Can anyone help me to understand why and how to avoid it ?

(1) sri is a NEDI function

%  2X enlarge

function y=sri(x,level)

     for l=1:level
     y1=sri1(x);y=sri2(y1);
     x=y;
     end

(2) sri1 is another NEDI function

function y=sri1(x)

[M,N]=size(x);
%e=zeros(M,N);
for i=1:2
   for j=1:2
      y(i:2:2*M,j:2:2*N)=x;
   end
end

T=5;cnt=1;
ix=[-1 -1 1 1];
iy=[-1 1 -1 1];
mx=[0 0 1 1];
my=[0 1 0 1];
cnt=1;
for i=-2:4
   for j=-2:4
      nx(cnt)=i;ny(cnt)=j;cnt=cnt+1;
  end
end
%nx=[0 1 1 0 -1 -1  0  1 2 2 0 1 -1 -1 2 2];
%ny=[0 0 1 1  0  1 -1 -1 0 1 2 2 -1  2 -1 2];

th=8;
for i=T:M-T
   for j=T:N-T
      for k=1:4
      C(:,k)=diag(x(i+ix(k)+nx,j+iy(k)+ny));
     end
      r=diag(x(i+nx,j+ny));
      s=diag(x(i+mx,j+my));
      if det(C'*C)==0|var(s)<th
         a=ones(4,1)/4;
      else
         a=inv(C'*C)*(C'*r);
         %a=fun(C,r);
      end
      y(2*i,2*j)=sum(a.*s); 
      if y(2*i,2*j)<0|y(2*i,2*j)>255
         a=ones(4,1)/4;y(2*i,2*j)=sum(a.*s);
    end
   end
end

(3) sri2 is also another NEDI function

% second step of edge-directed interpolation
function y=sri2(x)
y=x;
[M,N]=size(x);
T=8;cnt=1;
mx=[0 0 1 -1];
my=[-1 1 0 0];
%nx=[-1 0 0 1 -2 1 -1 2 1 2 -2 -1 3 -3 0 0]; 
%ny=[0 -1 1 0 1 -2 2 -1 2 1 -1 -2 0 0 -3 3];
cnt=1;
for i=-2:3
   for j=-2:3
      nx(cnt)=i+j-1;ny(cnt)=i-j;cnt=cnt+1;
   end
end

th=8;
for i=T:2:M-T
   for j=T+1:2:N-T
       for k=1:4
           C(:,k)=diag(x(i+2*mx(k)+nx,j+2*my(k)+ny));
       end
      r=diag(x(i+nx,j+ny));
      s=diag(x(i+mx,j+my));
      if det(C'*C)==0|var(s)<th
         a=ones(4,1)/4;
     else
          a=inv(C'*C)*(C'*r);
%          %a=fun(C,r);
end
      y(i,j)=sum(a.*s);
      if y(i,j)<0|y(i,j)>255
         a=ones(4,1)/4;y(i,j)=sum(a.*s);
    end
   end
end

for i=T+1:2:M-T
   for j=T:2:N-T
       for k=1:4
       C(:,k)=diag(x(i+2*mx(k)+nx,j+2*my(k)+ny));
   end
       r=diag(x(i+nx,j+ny));
      s=diag(x(i+mx,j+my));
      if det(C'*C)==0|var(s)<th
         a=ones(4,1)/4;
     else
         a=inv(C'*C)*(C'*r);
%          %a=fun(C,r);
end
      y(i,j)=sum(a.*s);
      if y(i,j)<0|y(i,j)>255
         a=ones(4,1)/4;y(i,j)=sum(a.*s);
    end
   end
end