function y=Mult_Dirac_W(x);
% Multiplies a vector by the Wilson-Dirac matrix
global N kp1 kp2 kp3 kp4 km1 km2 km3 km4
global g1_1 g2_1 g3_1 g4_1 g5_1 g6_1 g7_1 g8_1
global g1_2 g2_2 g3_2 g4_2 g5_2 g6_2 g7_2 g8_2
global g1_3 g2_3 g3_3 g4_3 g5_3 g6_3 g7_3 g8_3

mass=-.82;

t0=clock();
x=reshape(x,3*N,4);
x1=reshape(x(:,1),3,N);
x2=reshape(x(:,2),3,N);
x3=reshape(x(:,3),3,N);
x4=reshape(x(:,4),3,N);

% Dirac_Wilson_12
      y1_nn1 = x1(:,kp1) + i*x4(:,kp1);
      y2_nn1 = x2(:,kp1) + i*x3(:,kp1);
      y1_nn2 = x1(:,kp2) +   x4(:,kp2);
      y2_nn2 = x2(:,kp2) -   x3(:,kp2);
      y1_nn3 = x1(:,kp3) + i*x3(:,kp3);
      y2_nn3 = x2(:,kp3) - i*x4(:,kp3);
      y1_nn4 = x1(:,kp4) +   x3(:,kp4);
      y2_nn4 = x2(:,kp4) +   x4(:,kp4);

      y1_nn5 = x1(:,km1) - i*x4(:,km1);
      y2_nn5 = x2(:,km1) - i*x3(:,km1);
      y1_nn6 = x1(:,km2) -   x4(:,km2);
      y2_nn6 = x2(:,km2) +   x3(:,km2);
      y1_nn7 = x1(:,km3) - i*x3(:,km3);
      y2_nn7 = x2(:,km3) + i*x4(:,km3);
      y1_nn8 = x1(:,km4) -   x3(:,km4);
      y2_nn8 = x2(:,km4) -   x4(:,km4);
   
% Gauge field multiplication
y1_nn1=[sum(y1_nn1.*g1_1); sum(y1_nn1.*g1_2); sum(y1_nn1.*g1_3)];
y1_nn2=[sum(y1_nn2.*g2_1); sum(y1_nn2.*g2_2); sum(y1_nn2.*g2_3)];
y1_nn3=[sum(y1_nn3.*g3_1); sum(y1_nn3.*g3_2); sum(y1_nn3.*g3_3)];
y1_nn4=[sum(y1_nn4.*g4_1); sum(y1_nn4.*g4_2); sum(y1_nn4.*g4_3)];
y1_nn5=[sum(y1_nn5.*g5_1); sum(y1_nn5.*g5_2); sum(y1_nn5.*g5_3)];
y1_nn6=[sum(y1_nn6.*g6_1); sum(y1_nn6.*g6_2); sum(y1_nn6.*g6_3)];
y1_nn7=[sum(y1_nn7.*g7_1); sum(y1_nn7.*g7_2); sum(y1_nn7.*g7_3)];
y1_nn8=[sum(y1_nn8.*g8_1); sum(y1_nn8.*g8_2); sum(y1_nn8.*g8_3)];

y2_nn1=[sum(y2_nn1.*g1_1); sum(y2_nn1.*g1_2); sum(y2_nn1.*g1_3)];
y2_nn2=[sum(y2_nn2.*g2_1); sum(y2_nn2.*g2_2); sum(y2_nn2.*g2_3)];
y2_nn3=[sum(y2_nn3.*g3_1); sum(y2_nn3.*g3_2); sum(y2_nn3.*g3_3)];
y2_nn4=[sum(y2_nn4.*g4_1); sum(y2_nn4.*g4_2); sum(y2_nn4.*g4_3)];
y2_nn5=[sum(y2_nn5.*g5_1); sum(y2_nn5.*g5_2); sum(y2_nn5.*g5_3)];
y2_nn6=[sum(y2_nn6.*g6_1); sum(y2_nn6.*g6_2); sum(y2_nn6.*g6_3)];
y2_nn7=[sum(y2_nn7.*g7_1); sum(y2_nn7.*g7_2); sum(y2_nn7.*g7_3)];
y2_nn8=[sum(y2_nn8.*g8_1); sum(y2_nn8.*g8_2); sum(y2_nn8.*g8_3)];


% Dirac_Wilson_1234
% y = sum over mu=1,4 : [(1-gamma_mu)*y_kp + (1+gamma_mu)*y_km]
y1 =    y1_nn1 + y1_nn2 +   y1_nn3 + y1_nn4 +   y1_nn5 + y1_nn6 +   y1_nn7 + y1_nn8;
y2 =    y2_nn1 + y2_nn2 +   y2_nn3 + y2_nn4 +   y2_nn5 + y2_nn6 +   y2_nn7 + y2_nn8;
y3 = -i*y2_nn1 - y2_nn2 - i*y1_nn3 + y1_nn4 + i*y2_nn5 + y2_nn6 + i*y1_nn7 - y1_nn8;
y4 = -i*y1_nn1 + y1_nn2 + i*y2_nn3 + y2_nn4 + i*y1_nn5 - y1_nn6 - i*y2_nn7 - y2_nn8;

y=[y1(:);y2(:);y3(:);y4(:)];
y=(mass+4)*x(:)-0.5*y;
% Copyright (C) 2006-2007 Artan Borici.
% This program is a free software licensed under the terms of the GNU General Public License