To give a bit more detail on the calculations using Mathematica, this is an example of code which I used to calculate Duijvestijn's 112x112 Perfect Squared Square with 21 squares (the Perfect Squared Square with the fewest squares);
I have modified the code slightly to insert the graph of the GreatRhombicosidodecahedron in place of Duijvestijn's graph;
It appears the original error has disappeared and this code now works;
Text[" Incidence matrix Aa[i,j] -> graph G(i,j); i nodes, j edges
Aa[i,j]={ 1, if edge j is directed towards node i,
Aa[i,j]={-1, if edge j is directed from node i,
Aa[i,j]={ 0, if edge j is not incident to node i. "]
Text[" The Mathematica Archimedean Graphs are undirected and need to turned into directed graphs (any direction will work)]
Ga = DirectedGraph[GraphData["GreatRhombicosidodecahedralGraph"], "Random"]
Aa = IncidenceMatrix[Ga]
Text[" A row and column of the incidence matrix need to be removed"]
A = Drop[Aa, {1, 1}] ;
MatrixForm[A]
Text["Multiply A by its transpose, Transpose[A] to get the Laplacian
(or Kirchhoff) matrix K"]
K = A.Transpose[A];
MatrixForm[K]
Text["Perform LU decomposition on K, to factor it into an upper (u)
and a lower (l) triangular matrix and diagonal matrix"]
{lu, p, c} = LUDecomposition[K];
det = Tr[lu, Times]
Text["The V matrix is equal to determinant* Inverse matrix of K and
gives the (full)voltages for nodes satisfying Kirchhoff's 1st law"]
V = det*Inverse[K];
MatrixForm[V]
Text["Dividing rows or columns by GCD , we get reduced voltages -
this step is useful for calculating the Bouwkampcode "]
W = Apply[GCD, V, {1}]
Y = MatrixForm[V/W]
Text["The triple matrix product A.V.A gives the edge full current
solutions matrix F"]
F = Transpose[A].V.A;
MatrixForm[F]
Text["Dividing rows or columns of F by the GCD of each row we get a
reduction vector, needed to reduce the full currents so we get
tilings in smallest solutions"]
R = Apply[GCD, F, {1}]
Text["Dividing the full currents matrix F by reduction vector R gives
the reduced currents matrix B, which provides all the solutions in
terms of the sizes of the edge currents. Edge currents become square
sizes in the squared square. The last row of the B matrix is
Duijvestijn's 112 squared square, with all the square sizes in
reduced (GCD=1) size. The other rows are other squared
rectangles, and the non-diagonal entries are element (square) sizes, and the width of the rectangle is the
diagonal entry in each row, once we know the width we can work out
the height from the determinant. In the case of a squared square w =
h = diagonal entry."]
Text["Now we substitute the values of the edge currents from matrix
B, one row at a time, into the edges of the graph we started with,
this gives us the 'Smith Diagram' of the squared square, and we can
now draw or construct the squared square. A separate program to produce Bouwkampcode from the reduced currents and voltages is how this is usually done."]
B = MatrixForm[F/R]
