I wouldn't call it GL, it is the orthogonal group of the lattice we are discussing. References, as i said, include *Lattices and Codes* by W. Ebeling, *Rational Quadratic Forms* by Cassels, these two being available and inexpensive.  

We double the quadratic form to get an integral Gram matrix:

$$ G \; = \;  
 \left(  \begin{array}{rrr}
  2 & 1 & 0 \\\
  1  & 2 & 0 \\\
  0 & 0  & -2  
\end{array} 
  \right).
  $$

Your solutions, the light cone, are column vectors $v$ such that $ v^T G v = 0.$ A root, since we have doubled everything to get an "even" lattice, is a vector $r$ with norm 2, $ r^T G r = 2.$ The general definition of reflection (Cassels calls this a symmetry, page 19) in any vector $w$  is that 
$$ x \mapsto \; \; x \; -  \; \frac{2 \, x^T G w}{w^T G w} \; w.$$
As a result, when we take $w$ to be a root, the factors of 2 cancel and we are taking lattice points to other lattice points. In your original form, a root $(x,y,z)$ solves $x^2 + x y + y^2 = 1 + z^2.$ Then the reflection in the root is just a linear map, determinant $-1,$ and is therefore given by a square matrix with respect to the original basis. Finally, the reflection is an isometry, part of the orthogonal group of the quadratic form, and if we call the matix $A,$ it solves $$A^T G A = G.$$ 


I think it wise to include reflection in the root  $(1,0,0)^T$ to get some negative values taken care of,
$$ A_0 \; = \;  
 \left(  \begin{array}{rrr}
-1 & -1 & 0 \\\
0 & 1 & 0 \\\
0 & 0 & 1
\end{array} 
  \right).
  $$ 

We also include reflection in the root  $(0,1,0)^T,$ 
$$ A_{00} \; = \;  
 \left(  \begin{array}{rrr}
1 & 0 & 0 \\\
-1 & -1 & 0 \\\
0 & 0 & 1
\end{array} 
  \right).
  $$ 


Note that we have already included the automorph that interchanges the first two items, $$ (x,y,z) \mapsto (y,x,z). $$

$$ A_0 A_{00} A_0 \; = \;  
 \left(  \begin{array}{rrr}
0 & 1 & 0 \\\
1 & 0 & 0 \\\
0 & 0 & 1
\end{array} 
  \right).
  $$ 

Although it may not be called a root, having norm $-2,$ we include reflection in  $(0,0,1)^T,$ or
$$ A_{000} \; = \;  
 \left(  \begin{array}{rrr}
1 & 0 & 0 \\\
0 & 1 & 0 \\\
0 & 0 & -1
\end{array} 
  \right).
  $$ 

Also, as is usual, we explicitly include  negation,
$$ (x,y,z) \mapsto (-x,-y,-z), $$ so 
$$ A_{0000} \; = \;  
 \left(  \begin{array}{rrr}
-1 & 0 & 0 \\\
0 & -1 & 0 \\\
0 & 0 & -1
\end{array} 
  \right).
  $$ 

It occurred to me that we would get some reflections with quite small matrix entries by taking other vectors of norm $-2,$ for example $(1,1,2)^T.$
$$ A_{-1} \; = \;  
 \left(  \begin{array}{rrr}
4 & 3 & -4 \\\
3 & 4 & -4 \\\
6 & 6 & -7
\end{array} 
  \right).
  $$ This still solves, as with all the $A'$s, both $A^T G A = G$ and $A^2 = I.$ 

Taking the root $(3,4,6)^T$ gives the reflection 
$$ A_1 \; = \;  
 \left(  \begin{array}{rrr}
 -29 & -33 & 36 \\\
-40 & -43 & 48 \\\
-60 & -66 & 73 
\end{array} 
  \right).
  $$ 

If we take the initial triple in the light cone to be $c = (3,5,7)^T,$ we get $A_1 c = (0,1,1)^T$ which is indeed another solution, and shows that care must be used in constructing the "tree."



Taking the root $(5,15,18)^T$ gives the reflection 
$$ A_2 \; = \;  
 \left(  \begin{array}{rrr}
-124 & -175 & 180 \\\
-375 & -524 & 540 \\\
-450 & -630 & 649
\end{array} 
  \right).
  $$ 
we get $A_2 c = (13,35,43)^T$ 

 



Taking the root $(8,19,24)^T$ gives the reflection 
$$ A_3 \; = \;  
 \left(  \begin{array}{rrr}
-279 & -368 & 384 \\\
-665 & -873 & 912 \\\
-840 & -1104 & 1153
\end{array} 
  \right).
  $$ 
we get $A_3 c = (11,24,31)^T$ 




Taking the root $(7,32,36)^T$ gives the reflection 
$$ A_4 \; = \;  
 \left(  \begin{array}{rrr}
-321 & -497 & 504 \\\
-1472 & -2271 & 2304 \\\
-1656 & -2556 & 2593
\end{array} 
  \right).
  $$ 
we get $A_4 c = (80,357,403)^T$ 




Taking the root $(12,47,54)^T$ gives the reflection 
$$ A_5 \; = \;  
 \left(  \begin{array}{rrr}
-851 & -1272 & 1296 \\\
-3337 & -4981 & 5076 \\\
-3834 & -5724 & 5833
\end{array} 
  \right).
  $$ 
we get $A_5 c = (159,616,709)^T$ 

One typically includes $\pm 1$ anyway. Ian Agol would know how many reflections are enough, but I suspect this will do. Nothing really wrong with finding too many.
Note that the squarefree parts in $1+z^2$ in the five nontrivial roots I chose are $37, 13, 577, 1297, 2917.$