function A=Dirac4(u1,u2,u3,u4);
% Constructs Wilson-Dirac operator
mass=0; N1=8; N2=8; N3=8; N4=16; N=N1*N2*N3*N4;
% Gamma matrices
gamma1 = [0, 0, 0,-i;  0, 0,-i, 0;  0, i, 0, 0;  i, 0, 0, 0];
gamma2 = [0, 0, 0,-1;  0, 0, 1, 0;  0, 1, 0, 0; -1, 0, 0, 0];
gamma3 = [0, 0,-i, 0;  0, 0, 0, i;  i, 0, 0, 0;  0,-i, 0, 0];
gamma4 = [0, 0,-1, 0;  0, 0, 0,-1; -1, 0, 0, 0;  0,-1, 0, 0];
% Projection operators
P1_plus = eye(4)+gamma1; P1_minus=eye(4)-gamma1;
P2_plus = eye(4)+gamma2; P2_minus=eye(4)-gamma2;
P3_plus = eye(4)+gamma3; P3_minus=eye(4)-gamma3;
P4_plus = eye(4)+gamma4; P4_minus=eye(4)-gamma4;
% Shift operators
p1=[N1,1:N1-1]; p2=[N2,1:N2-1]; p3=[N3,1:N3-1]; p4=[N4,1:N4-1];
I1=speye(N1); I2=speye(N2); I3=speye(N3); I4=speye(N4);
T1=I1(:,p1); T2=I2(:,p2); T3=I3(:,p3); T4=I4(:,p4);
E1=spkron(I4,spkron(I3,spkron(I2,spkron(T1,speye(3)))));
E2=spkron(I4,spkron(I3,spkron(T2,spkron(I1,speye(3)))));
E3=spkron(I4,spkron(T3,spkron(I2,spkron(I1,speye(3)))));
E4=spkron(T4,spkron(I3,spkron(I2,spkron(I1,speye(3)))));
% Gauge Field configuration {u1, u2, u3, u4}: 9*N by 2 matrices
I_N=speye(N);
[I,J]=spfind(spkron(I_N,ones(3)));
U1=sparse(I,J,u1(:,1)+i*u1(:,2),3*N,3*N);
U2=sparse(I,J,u2(:,1)+i*u2(:,2),3*N,3*N);
U3=sparse(I,J,u3(:,1)+i*u3(:,2),3*N,3*N);
U4=sparse(I,J,u4(:,1)+i*u4(:,2),3*N,3*N);
% Upper triangular
U=spkron(P1_minus,U1*E1)+spkron(P2_minus,U2*E2)+spkron(P3_minus,U3*E3)+spkron(P4_minus,U4*E4);
% Lower triangular
L=spkron(P1_plus ,U1*E1)+spkron(P2_plus ,U2*E2)+spkron(P3_plus ,U3*E3)+spkron(P4_plus ,U4*E4);
%M=U+L';
A=(mass+4)*speye(12*N)-0.5*(U+L');
% Copyright (C) 2006-2007 Artan Borici.
% This program is a free software licensed under the terms of the GNU General Public License