Expressions enabling systematic compilation of Hamiltonian and overlap matrix elements for an antisymmetrized multiterm geminal product trial function are derived, using double coset (DC) decompositions and subgroup adapted irreducible representations of the symmetric group, SN. The trial function may describe an even electron atomic or molecular system in any total spin eigenstate, and the geminals may be nonorthogonal, have arbitrary permutational symmetry, and be explicit functions of interelectronic distance. A DC decomposition is used to factor out permutations not exchanging particle labels between geminals (elements of the interior pair group, Sn2, n≡N/2). This reduces the sum over N! permutations to a sum over DC generators. If the irreducible representation λ (S) of SN is adapted to Sn2 each geminal is projected into its singlet or triplet component. The DC generators are chosen such that each has the form QP, where Q permutes odd particle labels only and P is a permutation of geminals (element of the exterior pair group, S?). With the aid of matrices called DC symbols an algorithm for these generators is derived, and used to find explicit sets for N=2, 4, 6, and 8. The N-electron Hamiltonian and overlap integrals arising with a particular DC generator QP are factored into products of smaller integrals, called cluster integrals, according to the cycle structure of Q. The cluster integrals are of only three main types—overlap, one-cluster energy, and two-cluster energy (analogs of orbital overlap, 1-electron, and 2-electron integrals) —and are further classified by order (number of geminals), geminal permutational symmetry, and in some cases pattern of connection. Matrix element compilation is systematic in that all N-electron integrals are products of a relatively small number of different types of cluster integral, and that N-electron integrals with similar factored forms are collected together in the summation. From counting the different cluster integrals required, it is concluded that a geminal product calculation not using orbital expansion is feasible only for systems with eight or less electrons. In some cases semiempirical calculations with correlated geminals might be considered, for the more complex cluster integrals (those of high order) are quite small for a system approximating a collection of localized electron pairs. The matrix element expressions are specialized for three cases—all geminals being singlets, strongly orthogonal geminals, and identical geminals. Comparison is made with a recently developed diagrammatical method.