Matrix algebra over finite fields can be implemented using the following procedure,
GFM:=proc(G)
local x,p,T,C;
global `&+`,`&*`,`&.`,`&^`,tr,det,ad,inv,Id,zero,
CharPoly,RandMat,RandVec;
uses LinearAlgebra;
x:=G:-variable;
p:=op([2,1,1],ifactors(G:-size));
T:=a->map(c->modp1(ConvertOut(c,x),p),a);
`&+`:=(a,b)->map(G:-ConvertIn,T(a)+T(b));
`&*`:=(a,b)->map(G:-`*`,b,a);
`&.`:=proc(a,b) C:=T(a).T(b);
if type(C,'rtable') then map(G:-ConvertIn,C)
else G:-ConvertIn(C) fi end;
tr:=G:-ConvertIn@Trace@T;
det:=G:-ConvertIn@Determinant@T;
ad:=c->map(G:-ConvertIn,Adjoint(T(c)));
inv:=c->G:-inverse(det(c)) &* ad(c);
Id:=()->Matrix(args,(i,j)->if i=j then G:-one else G:-zero fi);
zero:=()->Matrix(args,()->G:-zero);
`&^`:=(a,n::integer)->if n>0 then map(G:-ConvertIn,T(a)^n) elif n<0
then map(G:-ConvertIn,T(inv(a))^(-n)) else Id(op(1,a)) fi;
CharPoly:=(a,t)->map(c->modp1(Rem(c,G:-extension),p),
modp2(ConvertIn(CharacteristicPolynomial(T(a),t),t,x),p));
RandMat:=()->Matrix(args,G:-random);
RandVec:=n->Vector[`if`(type(procname,'indexed'),
op(procname),'column')](n,G:-random);
NULL end:

For example, G:=GF(3,4):
GFM(G);
M:=RandMat(3);
[ 2 3 3 2 ]
[1 + T + 2 T 1 + T 2 + 2 T ]
[ ]
M := [ 2 3 ]
[ 1 T + T + T 2 + T ]
[ ]
[ 2 3 2 3]
[ 1 + T + T 1 + 2 T + T 1 + 2 T + T + T ]
a:=G:-random();
2 3
a := 1 + T + T + T
a &* M;
[ 2 2 3 2 3]
[2 T + 2 T , 1 + T + 2 T , 2 T + T + 2 T ]
[ ]
[ 2 3 2 3 2]
[1 + T + T + T , 1 + T + T + 2 T , 1 + T + 2 T ]
[ ]
[ 2 3 2 3 2 3]
[2 + 2 T + 2 T , 2 + T + T + T , 1 + T + 2 T + T ]
M &+ M;
[ 2 3 3 2]
[2 + 2 T + T , 2 + 2 T , 1 + T ]
[ ]
[ 2 3 ]
[2 , 2 T + 2 T + 2 T , 1 + 2 T]
[ ]
[ 2 3 2 3]
[2 + 2 T + 2 T , 2 + T + 2 T , 2 + T + 2 T + 2 T ]
ad(M);
[ 3 2 3 ]
[2 + 2 T + T 2 + T + 2 T 0 ]
[ ]
[ 2 2 3 2 3]
[1 + T + 2 T 1 + T + T + 2 T 2 + 2 T + T ]
[ ]
[ 2 3 2 ]
[2 T + 2 T 1 + T 1 + 2 T + T ]
inv(M);
[ 2 3 ]
[2 T + 2 T 2 + 2 T 0 ]
[ ]
[ 3 3 2 3 ]
[ T + 2 T 2 + 2 T + 2 T 2 + T + T + 2 T ]
[ ]
[ 2 2 3]
[ 2 + T 2 + 2 T 2 + 2 T + 2 T + 2 T ]
inv(M) &. M;
[1 0 0]
[ ]
[0 1 0]
[ ]
[0 0 1]
CharPoly(M,t);
3 2 3 2 3
(1 + 2 T ) + (1 + 2 T + 2 T ) t + (1 + 2 T ) t + (1) t
det(M);
3
2 + T
tr(M);
3
2 + T
Id(3);
[1 0 0]
[ ]
[0 1 0]
[ ]
[0 0 1]
Id(3,5);
[1 0 0 0 0]
[ ]
[0 1 0 0 0]
[ ]
[0 0 1 0 0]
zero(2);
[0 0]
[ ]
[0 0]
zero(2,3);
[0 0 0]
[ ]
[0 0 0]
M &^ 2;
[ 2 3 2 ]
[ 2 + T + T 2 + 2 T + T 1 + T ]
[ ]
[ 2 3 2 3 2]
[ T + 2 T + T T + T 2 + T ]
[ ]
[ 2 3 2 3 ]
[1 + 2 T + T + 2 T 2 + T + T 1 + T ]
M &^ (-2);
[ 2 ]
[ 2 + T + 2 T 2 + T T ]
[ ]
[ 2 2 3 ]
[ 2 + T + 2 T 1 + 2 T 2 + T + T ]
[ ]
[ 2 3 3 2 3]
[2 T + T + 2 T 2 + T + 2 T 1 + 2 T + T + 2 T ]
M &^ 0;
[1 0 0]
[ ]
[0 1 0]
[ ]
[0 0 1]
b:=RandVec(3);
[ 2]
[T + 2 T ]
b := [ ]
[ 2 + T ]
[ ]
[ 1 ]
M &. b;
[ 3 ]
[ 2 + 2 T ]
[ ]
[ 2]
[1 + 2 T + T ]
[ ]
[ 2 T ]
c:=RandVec[row](3);
[ 2 3 3 2]
c := [2 + T + T , 2 + T + 2 T , 2 + 2 T + 2 T ]
c &. M;
[ 3 2 3 2 3]
[2 + T , 2 + 2 T + T + T , 2 + 2 T + 2 T + T ]
c &. b;
2
1 + 2 T + T
A:=b &. c;
[ 3 2 3 2 ]
[ 2 T 1 + T + T 2 + 2 T ]
[ ]
A := [ 2 2 3 3 ]
[ T 2 + 2 T + T 1 + 2 T ]
[ ]
[ 2 3 3 2]
[2 + T + T 2 + T + 2 T 2 + 2 T + 2 T ]

Later I will add a few more commands. The purpose of this post is to provide functionality (in a short and simple way), not effectiveness. Effective procedures should be written in assembler anyway, not in Maple. An example of that can be seen in Binomial coefficients mod 2.