The luminosity function of galaxy systems - from single galaxies and small groups to rich clusters - is determined. I first determine the rich-cluster luminosity function, then a “groupings” luminosity function, which includes all systems from small groups to rich clusters, and finally, I estimate a general luminosity function of all galaxy systems (AGS), from single galaxies to rich clusters. The catalogs of Abell (rich clusters) and Turner and Gott (groups, singles) are used. The AGS luminosity function is found to be a smooth and steeply decreasing function of number density with luminosity; it spans over twelve decades of number densities in the observed five-decade range of luminosities (~109 to 1014 L⊙). The function exhibits an exponential cutoff in the rich-cluster luminosity domain, a power law dependence over most of the intermediate luminosity range, and a flattening of slope below the characteristic galactic luminosity, L*. A smooth transition between the groups and rich-clusters occurs near ~1013L⊙. The groupings LF (i.e., groups to rich-clusters) can be represented well by a Schechter-type analytic form η(L)=ηo (L/Lo)−α exp (-L/Lo), where η(L) is number of systems per unit volume per unit luminosity. The best-fit parameters are α=2.0 ± 0.1, Lo=1.0 ± 0.2×1013 L⊙, and a normalization ηo=1.6 × 10−7 systems Mpc−3 (1012L⊙)−1. The “break” at Lo≃1 × 1013L⊙ occurs at the typical luminosity of rich Abell-type clusters. The slope α is found to depend on the density enhancement factor with which the groups are selected, but the dependence is relatively weak. The exponential decay at the bright end (rich clusters) has an equivalent power-law slope of approximately −5, i.e., η(L)∝L−5. The functional form of the LF is consistent with the Press-Schechter (1974) model of galaxy formation, and with the recent Silk-White (1978) model (predicted slope of −2.0). When single galaxies are added to the groupings luminosity function, the luminosity function of all galaxy systems (AGS) can be roughly described as η(L)=1.6 × 10−7 (L/Lo)−2e-L/Lo (1 + 1.6 × 1010/L)−0.75Mpc−3(1012L⊙)−1 with Lo=1 × 1013L⊙. This function approaches Schechter's field galaxy luminosity function at low luminosities, and the groupings function at high luminosities. Comparisons of the LF with the LF of galaxies and X-ray clusters of galaxies are discussed, and an estimate of is given for the observed luminosity density in the universe.The above luminosity function of galaxy systems could be compared with high redshift cluster data, when available, and thus shed light on possible evolutionary processes.The X-ray luminosity function of clusters of galaxies is calculated using the general LF of optical clusters above and a relation between X-ray and optical luminosity based on a thermal bremmstrahlung model. The derived X-ray LF applies to clusters whose X-ray emission is mostly due to thermal bremmstrahlung from a smooth, hot intracluster gas. The main features predicted by the LF are a change of slope in the log ηx -log Lx plane at Lx ~1044 ergs/s, a moderately steep slope of −2.5 at high luminosities (~1044-1046 ergs/s), and a flatter slope of ~-1.3 below 1044 ergs/s (down to ~1042 or 1043 ergs/s, where the smooth thermal bremsstrahlung may not dominate the total X-ray emission). Simple analytic expressions of the LF are given for various X-ray energy bands. X-ray data from Uhuru and Ariel-5, available for the 1044-1046 ergs/s luminosity range (including upper limits in the 1045-1046 range), agree well with the LF calculated here. The present calculations suggest that at least to Lx ~ 1046 ergs/s the X-ray LF does not decay faster than ~Lx−2.5 (for a wide enough energy band). It is also suggested that fewer faint X-ray clusters than previously expected may exist (and be detected by the Einstein satellite).