Permutability Graph with GAP Could someone suggest a piece of code that will produce
the subgroup permutability graph of a finite group? There is more than one way this
graph can be defined. The one I have in mind has the non-normal subgroups of the group as vertices and two vertices are joined by an edge if the corresponding subgroups permute.
 A: First, let's pick up some group, for example this one, which is a semidirect product of $C_9 \times C_3$ and $C_3$:
gap> G:=SmallGroup(81,9);
<pc group of size 81 with 4 generators>
gap> StructureDescription(G);
"(C9 x C3) : C3"

The next command calculates the list of all its non-normal subgroups by iterating over each conjugacy class of subgroups, containing more than 1 subgroup:
gap> cc:=Concatenation( List( Filtered( ConjugacyClassesSubgroups(G),
>                             c -> Size(c)>1), AsList ) )
[ Group([ f3 ]), Group([ f3*f4 ]), Group([ f3*f4^2 ]), Group([ f1 ]), ...
... some more lines omitted ...
... Group([ f1*f2^2*f3^2, f4 ]), Group([ f1*f2^2*f3, f4 ]) ]
gap> Length(cc);
42

Remarkably, it has 42 of them, what makes this group an excellent example :-) 
The next function constructs the adjacency matrix for the permutability graph. It first creates the identity matrix (obviously, each subgroup is permutable with itself, so we have $1$s on the main diagonal and do not have to check this), and then populates the rest of the entries:
PermutabilityMatrix:=function(cc)
local A, i, j;
A := IdentityMat(Length(cc));
for i in [1..Length(cc)-1] do
  for j in [i+1..Length(cc)] do 
    if Size(ClosureGroup(cc[i],cc[j]))*Size(Intersection(cc[i],cc[j])) = 
       Size(cc[i])*Size(cc[j]) then 
      A[i][j]:=1;
      A[j][i]:=1;
    fi;
  od; 
od;
return A;
end;

One could do even better with the Permut package:
LoadPackage("permut");
PermutabilityMatrixWithPERMUT:=function(cc)
local A, i, j;
A := IdentityMat(Length(cc));
for i in [1..Length(cc)-1] do
  for j in [i+1..Length(cc)] do 
    if ArePermutableSubgroups(cc[i],cc[j]) then
      A[i][j]:=1;
      A[j][i]:=1;
    fi;
  od; 
od;
return A;
end;

which have some shortcuts for partial cases. Just to check that both functions agree about the result:
gap> A:=PermutabilityMatrix(cc);;
gap> B:=PermutabilityMatrixWithPERMUT(cc);;
gap> A=B;
true

Now you may load the GRAPE package and construct GRAPE graph given by this adjacency matrix with the commmands:
LoadPackage("grape");
g:=Graph( Group(()), [1..Length(cc)], OnPoints, function(x,y) return A[x][y]=1;end,true);

This allows you to explore properties of this graph, for example:
gap> Diameter(g);
4
gap> IsConnectedGraph(g);
true

One could further call AssignVertexNames(g,cc); to label vertices of the graph with subgroups.

Update:. Trying to understand and reproduce example from page 6 of the paper you've mentioned, I've made some quickly-written code, which may be useful to you for further experiments. First, this function returns a list of all non-normal subgroups:
AllNonNormalSubgroups:= function(G)
local c;
return Concatenation( List( Filtered( ConjugacyClassesSubgroups(G),
                            c -> Size(c)>1), AsList ) );
end;

Remember that their ordering is not guaranteed to be the same for different instances of $G$ as probabilistic methods may be involved.
Next, GAP operates with binary relations - for example, finds a transitive closure. This function looks similar to adjacency matrix, but it returns an input for BinaryRelationOnPoints:
PermutabilityRelation:=function(cc)
local r, i, j;
r := List( [ 1 .. Length(cc) ], i -> [i]);
for i in [1..Length(cc)-1] do
  for j in [i+1..Length(cc)] do 
    if Size(ClosureGroup(cc[i],cc[j]))*Size(Intersection(cc[i],cc[j])) = 
       Size(cc[i])*Size(cc[j]) then 
      AddSet(r[i],j);
      AddSet(r[j],i);
    fi;
  od; 
od;
return r;
end;

In the following example we will show that $C_3 \times S_3$ is irreducible with respect to permutability (as defined in the paper) since the transitive closure of the permutability relation has one equivalence class:
gap> G:=DirectProduct(CyclicGroup(3),SymmetricGroup(3));
gap> nns := AllNonNormalSubgroups(G);;
gap> pr := PermutabilityRelation(nns);
[ [ 1, 6 ], [ 2, 7 ], [ 3, 8 ], [ 4, 5, 6, 7, 8 ], [ 4, 5, 6, 7, 8 ], 
  [ 1, 4, 5, 6 ], [ 2, 4, 5, 7 ], [ 3, 4, 5, 8 ] ]
gap> R  := BinaryRelationOnPoints( pr );
Binary Relation on 8 points
gap> T  := TransitiveClosureBinaryRelation ( R );
Binary Relation on 8 points
gap> IsEquivalenceRelation( T );
true
gap> EquivalenceClasses( T );
[ {1} ]

The following function takes a binary relation and returns a GRAPE graph corresponding to it:
LoadPackage("grape");
PermutabilityGraphByRelation:=function( rel )
local n;
n := DegreeOfBinaryRelation(rel);
return Graph( Group(()), [1..n], OnPoints, 
       function(x,y) return y in Successors(rel)[x]; end, true );
end;

Now we may explore graphs given by R and T above:
gap> PermutabilityGraphByRelation( R );;
gap> Diameter(last);
4
gap> PermutabilityGraphByRelation( T );;
gap> Diameter(last);
1

On the other hand, for $D_{32}$ we have two classes:
gap> G:=DihedralGroup(32);
<pc group of size 32 with 5 generators>
gap> nns := AllNonNormalSubgroups(G);;
gap> pr := PermutabilityRelation(nns);;
gap> R  := BinaryRelationOnPoints( pr );;
gap> T  := TransitiveClosureBinaryRelation ( R );;
gap> EquivalenceClasses( T );
[ {1}, {9} ]

GRAPE can construct graph as above, but it can find connected components only if the graph is simple. So if we want to find the diameter of each of the two connected components, we have to construct two graphs first, and then find their diameters. Luckily, GRAPE already has functionality for that via InducedSubgraph:
gap> pg:=PermutabilityGraphByRelation(R);;
gap> erp:=EquivalenceRelationPartition(T);
[ [ 1, 2, 3, 4, 5, 6, 7, 8, 17, 18, 19, 20, 25, 26 ], 
  [ 9, 10, 11, 12, 13, 14, 15, 16, 21, 22, 23, 24, 27, 28 ] ]
gap> List(erp, w -> Diameter(InducedSubgraph(pg,w)));
[ 2, 2 ]

Thus, both of the connected components for permutability graph has the diameter 2, so $\Delta(G_{32})=2$.

Update (communicated to me by Leonard Soicher)
The code posted above may be improved to use the full power of GRAPE package to associate
a group acting as a group of automorphisms of the graph we are constructing, as illustrated below:
gap> LoadPackage("grape");
-----------------------------------------------------------------------------
Loading  GRAPE 4.6.1 (GRaph Algorithms using PErmutation groups)
by Leonard H. Soicher (http://www.maths.qmul.ac.uk/~leonard/).
Homepage: http://www.maths.qmul.ac.uk/~leonard/grape/
-----------------------------------------------------------------------------
true
gap> G:=SmallGroup(81,9);; # for example
gap> nonnormalsubgroups:=Concatenation( List( Filtered( ConjugacyClassesSubgroups(G),
>                                             c -> Size(c)>1), AsList ) );;
gap> permutable := function(A,B)
> return Size(ClosureGroup(A,B))*Size(Intersection(A,B)) = Size(A)*Size(B);
> end;
function( A, B ) ... end

Now nonnormalsubgroups is a $G$-invariant list, and permutable defines a $G$-invariant binary relation on nonnormalsubgroups (for the conjugation action of $G$ on its subgroups). Hence the required GRAPE graph can be computed as:
gap> gamma:=Graph(G,nonnormalsubgroups,OnPoints,permutable,true);;
gap> Diameter(gamma);
4

If you want this graph with the loops removed so as to obtain a simple graph, you can do:
gap> delta:=Graph(G,nonnormalsubgroups,OnPoints,
> function(A,B) return A<>B and permutable(A,B); end,true);;
gap> Diameter(delta);
4
gap> Girth(delta);
3

A: I think there can be problems if you define the graph with conjugacy classes of non-normal subgroups instead of considering all subgroups. For instance, take the semidirect product $G=[V]S$ of the symmetric group of degree $3$, $S=\langle a,b\mid a^3=b^2=1\rangle $ say, by a faithful and irreducible module $V=\langle c,d\rangle$ over the field of $7$ elements. This is
    SmallGroup(294, 7)

in GAP. Then $\langle b\rangle$ permutes with $\langle a\rangle$, but $\langle b\rangle$ does not permute with $\langle a^c\rangle$. Therefore I think that we will need all non-normal subgroups in the construction of the graph, as Alexander Konovalov says, not only the conjugacy classes of such subgroups.
