The setting is clear, but the question somehow interferes with the mathematical research question of finding a good $S$, subset in $\mathbb{F}_{p^2}$, $p$ prime.

My structural problem is that $S$ should be good for all entries in a matrix to be validated.
The property is defined elementwise, so we need the matrix with its entries.
(Some triangular or "parabolic" shapes are eliminated.)

Now to the questions:

(1) which matrix operations are available for elements of GU(n,q)?

I tried:

```
sage: G = GU( 3, 7 )
sage: G
General Unitary Group of degree 3 over Finite Field in a of size 7^2
sage: G.order().factor()
2^10 * 3 * 7^3 * 43
sage: A = G.random_element()
sage: A
[ a + 1 5 4*a + 1]
[ 6*a 6*a + 3 6*a]
[5*a + 6 2*a + 1 3*a + 1]
sage: A.
A.N A.conjugacy_class A.inverse A.matrix A.parent A.save
A.base_extend A.db A.is_idempotent A.multiplicative_order A.powers A.subs
A.base_ring A.dump A.is_one A.n A.pseudo_order A.substitute
A.cartesian_product A.dumps A.is_zero A.numerical_approx A.rename A.word_problem
A.category A.gap A.list A.order A.reset_name
```

Above, after A. i was hitting the TAB (twice), and the list of implemented (public) methods was printed.
An other way of investigating A is by printing dir(A).
The "hidden" methods show themselves after A._TABULATOR . (In the interpreter.)
The we can use for instance the method order.

```
sage: A.order().factor()
2^3 * 43
```

Yes, determinant is not in the landscape. I tried to coerce, to apply a morphism.
Finally, the following work around could indeed walk around.

```
p = 7
G = GU( 3, p )
F.<a> = GF( p^2 )
H = GL( 3, F )
H_Matrices = [] # and we will append.
for A in G.gens():
HA = H ( A.gap().sage() ) # tricky coercion, not proud of it
H_Matrices.append( HA )
HG = H.subgroup( H_Matrices )
print "HG has order: %s" % HG.order().factor()
```

And we get:

```
HG has order: 2^10 * 3 * 7^3 * 43
```

And in HG we also have the determinant. For instance:

```
h = HG.random_element()
print "The matrix\n%s\nhas determinant %s" % ( h, h.matrix().det() )
```

This gives a value in $\pm 1$:

```
The matrix
[ a 6*a + 1 4*a + 1]
[5*a + 5 3*a 2*a + 2]
[2*a + 3 3*a + 3 a + 3]
has determinant 6
```

(2)
For the second question, it is easy to implement a validator.
For instance:

```
A, B = H_Matrices
A = A.matrix()
B = B.matrix()
S = set( A.list() + B.list() )
# S is the list of entries in A, B
FA = H( [ F.frobenius_endomorphism()( entry ) for entry in A.list() ] ).matrix()
FB = H( [ F.frobenius_endomorphism()( entry ) for entry in B.list() ] ).matrix()
```

So we obtain the $F'$ action on matrices, which invokes the Frobenius $F$ and the reversion of rows and columns,
that invariates $A$.

Now we can build products, and check if the entries land in S.