Haskell, 79 bytes

Basically a port of ngenisis's Mathematica method. (Except I'm using 0-indexed arrays.)

c takes a list of lists of Ints and returns an integer.

c g=sum[1|l<-foldr(\r->([id,(r!!0:)]<*>))[[]]g,and[g!!x!!y`elem`l|x<-l,y<-l]]-1

Try it online!

It is assumed that the Ints are numbered the same as the rows (outer lists) and columns showing their multiplication. Thus, since 0 is the identity, the first column is the same as the indices of the rows. This allows using the entries of the first column to construct the subsets.

How it works

• c is the main function.
• g is the group array as a list of lists.
• l is a subset of the elements. The list of subsets is constructed as follows:
  • foldr(\r->([id,(r!!0:)]<*>))[[]]g folds a function over the rows of g.
  • r is a row of g, whose first (0th) element is extracted as an element that may be included ((r!!0:)) or not (id).
  • <*> combines the choices for this row with the following ones.
• and[g!!x!!y`elem`l|x<-l,y<-l] tests for each pair of elements in l whether their multiple is in l.
• sum[1|...]-1 counts the subsets that pass the test, except for one, the empty subset.