Adjacency matrix of a hexagonal graph

I am looking for an efficient way to compute the adjacency matrix of a hexagonal graph, so that away from the boundary the graph is 3-regular and its minimal cycles are hexagons. The graph is arranged in a rectangle with n vertices in each row and m rows. My follow on computations require me to represent this matrix as a sum of three separate adjacency matrices: one for the horizontal edges, one for the edges that slant in the 'North West' direction and one for the North East. I give my code below; I need to make the whole computation much quicker (this is one of many elements), so I would welcome any suggestions. Thank you.

%% Hexagonal graph
% This function produces the adjacency matrix of a "rectangular graph" whose cells are hexagons. 
% The size of the graph is [n,m], where n is the number of vertices
% on one level and m is the number of levels. The adjacency matrix is 
% given in three parts: 
% H: the horizontal edges, 
% and, looking from the bottom: 
% NE: edges slanted "up and right", 
% NW: edges slanted "up and left"

function [H,NW,NE]  = HexGraphPart(n,m)

GridH  = zeros(n*m);  
GridNW = zeros(n*m); 
GridNE = zeros(n*m); 

% NW part: 
% odd numbered rows: 
for kk = 1:floor(m/2)
    for ii = 1:2:n
        if (2*kk-1)*n + ii <= n*m
            GridNW((2*(kk-1))*n+ii, (2*kk-1)*n + ii) = 1; 
        end
    end
end
% even numbered rows: 
for kk = 1:floor(m/2)
    for ii = 2:2:n
        if 2*kk*n +ii  <= n*m
            GridNW((2*kk-1)*n + ii , 2*kk*n + ii ) = 1; 
        end
    end
end

% NE part:
% Odd numbered rows: 

for kk = 1:floor(m/2)
    for ii = 2:2:n
        if (2*kk-1)*n + ii <= n*m
            GridNE(2*(kk-1)*n+ii, (2*kk-1)*n + ii) = 1; 
        end
    end
end

% Even numbered rows: 

for kk = 1:floor(m/2)
    for ii = 1:2:n
        if 2*kk*n +ii  <= n*m
            GridNE((2*kk-1)*n + ii , 2*kk*n + ii ) = 1; 
        end
    end
end

% Horizontal edges: 
% Odd numbered rows
for kk = 1: ceil(m/2) 
    for ii = 1:2:n 
        if 2*(kk-1) +ii +1 <= n*m
        GridH(2*(kk-1)*n +ii, 2*(kk-1)*n +ii +1) = 1; 
        end
    end
end

% Even numbered rows: 

for kk = 1:floor(m/2)
    for ii = 2:2:n-1
        if (2*kk-1)*n + ii + 1 <= n*m
            GridH((2*kk-1)*n + ii, (2*kk-1)*n + ii + 1) = 1; 
        end
    end
end
    H =  GridH + transpose(GridH); 
    NW = GridNW + transpose(GridNW); 
     NE = GridNE + transpose(GridNE); 
end

Answer 1

[Bruno Luong]

[A,B,C] = HexGraphPart2(7,8);
G = graph(A+B+C); 
figure;
plot(G);

function [H,NW,NE]=HexGraphPart2(m,n)
[I,J]=ndgrid(1:m,1:n);
K = I + (J-1)*m;
p=m*n;
H = sparse(K(1:2:end-1,1:2:end),K(2:2:end,1:2:end),1,p,p) + ...
    sparse(K(2:2:end-1,2:2:end),K(3:2:end,2:2:end),1,p,p);
NW = sparse(K(1:2:end,1:2:end-1),K(1:2:end,2:2:end),1,p,p) + ...
     sparse(K(2:2:end,2:2:end-1),K(2:2:end,3:2:end),1,p,p);
NE = sparse(K(1:2:end,2:2:end-1),K(1:2:end,3:2:end),1,p,p) + ...
     sparse(K(2:2:end,1:2:end-1),K(2:2:end,2:2:end),1,p,p);
H=H+H';
NW=NW+NW';
NE=NE+NE';
end

NOTE:1 Your code is buggy with certain combination of m/n

NOTE2: the code can be simplied !if you only need A+B+C

[Jacek Brodzki]

Thank you for your suggestions. What combinations of m and n would create errors in my code? Regarding your Note 2, I do need separate adjacency matrices as I need to create operators from them in different ways, depending on the case.

[Bruno Luong]

Please try with (m,n)=(3,4) plot the graph, you'll see all the faces are pentagonal, rectangle, triangles, anything but hexagonal.

Correct graph

[Bruno Luong]

If one prefers to use "standard" hexagonal-pattern coordinates, and avoid graph-plotting method to decide on its own

m = 7;
n = 8;

[A,B,C,X,Y] = HexGraphPart2(m,n);
G = graph(A+B+C); 

close all
figure;
h = plot(G);
% Change node coordinates on regular hexagonal pattern
set(h,'XData',X(:),'YData',Y(:));
axis('equal');

function [H,NW,NE,X,Y]=HexGraphPart2(m,n)
x = (0:m-1)';
y = (0:n-1);
K = 1 + x + m*y;
p = m*n;
H  = sparse(K(1:2:end-1,1:2:end),K(2:2:end,1:2:end),1,p,p) + ...
     sparse(K(2:2:end-1,2:2:end),K(3:2:end,2:2:end),1,p,p);
NW = sparse(K(1:2:end,1:2:end-1),K(1:2:end,2:2:end),1,p,p) + ...
     sparse(K(2:2:end,2:2:end-1),K(2:2:end,3:2:end),1,p,p);
NE = sparse(K(2:2:end,1:2:end-1),K(2:2:end,2:2:end),1,p,p) + ...
     sparse(K(1:2:end,2:2:end-1),K(1:2:end,3:2:end),1,p,p);
H = H+H';
NW = NW+NW';
NE = NE+NE';
if nargout >= 5
    X = (3*x - mod(x+y,2))/2;
    Y = (sqrt(3)/2)*y + 0*x;
end
end

[Jacek Brodzki]

Thank you very much for this suggestion, Bruno. I have my own code to do compute the coordinates, but yours is much simpler. One question, why do you need

0*x

in the definition of Y?

[Bruno Luong]

Lower case y is a horizontal vector, x is vertical vector. Upper-case Y depends only y, but I need it as 2D array, thus + 0*x to expand it. You could use repmat, repelem as well if the trick perturbs you.

[Jacek Brodzki]

This is really neat; thank you again for taking the time to respond. This helped a lot.

Answer 2

[Jacek Brodzki]

Yes, you are right; I spotted that, too, but my main objective was to compute these objects for large values of m and n. Thank you for your suggestions, and taking the time to respond.

[Bruno Luong]

I believe the bug is related to odd/even characteristics of m and n, and not because they are small or large. I haven't try to look carefully your code.