The basic pam algorithm is fully described in chapter 2 of
  Kaufman and Rousseeuw(1990).  Compared to the k-means approach in kmeans, the
  function pam has the following features: (a) it also accepts a
  dissimilarity matrix; (b) it is more robust because it minimizes a sum
  of dissimilarities instead of a sum of squared euclidean distances;
  (c) it provides a novel graphical display, the silhouette plot (see
  plot.partition) (d) it allows to select the number of clusters
  using mean(silhouette(pr)[, "sil_width"]) on the result
  pr <- pam(..), or directly its component
  pr$silinfo$avg.width, see also pam.object.
When cluster.only is true, the result is simply a (possibly
  named) integer vector specifying the clustering, i.e.,
  pam(x,k, cluster.only=TRUE) is the same as 
  pam(x,k)$clustering but computed more efficiently.
The pam-algorithm is based on the search for k
  representative objects or medoids among the observations of the
  dataset.  These observations should represent the structure of the
  data.  After finding a set of k medoids, k clusters are
  constructed by assigning each observation to the nearest medoid.  The
  goal is to find k representative objects which minimize the sum
  of the dissimilarities of the observations to their closest
  representative object.
  
  By default, when medoids are not specified, the algorithm first
  looks for a good initial set of medoids (this is called the
  build phase).  Then it finds a local minimum for the
  objective function, that is, a solution such that there is no single
  switch of an observation with a medoid (i.e. a ‘swap’) that will
  decrease the objective (this is called the swap phase).
When the medoids are specified (or randomly generated), their order does not
  matter; in general, the algorithms have been designed to not depend on
  the order of the observations.
The pamonce option, new in cluster 1.14.2 (Jan. 2012), has been
  proposed by Matthias Studer, University of Geneva, based on the
  findings by Reynolds et al. (2006) and was extended by Erich Schubert,
  TU Dortmund, with the FastPAM optimizations.
The default FALSE (or integer 0) corresponds to the
  original “swap” algorithm, whereas pamonce = 1 (or
  TRUE), corresponds to the first proposal .... 
  and pamonce = 2 additionally implements the second proposal as
  well. 
The key ideas of ‘FastPAM’ (Schubert and Rousseeuw, 2019) are implemented
  except for the linear approximate build as follows:
    - pamonce = 3:
- reduces the runtime by a factor of O(k) by exploiting
      that points cannot be closest to all current medoids at the same time. 
    
- pamonce = 4:
- additionally allows executing multiple swaps
      per iteration, usually reducing the number of iterations. 
    
- pamonce = 5:
- adds minor optimizations copied from the
      - pamonce = 2approach, and is expected to be the fastest of the
      ‘FastPam’ variants included.
 
  
‘FasterPAM’ (Schubert and Rousseeuw, 2021) is implemented via
    - pamonce = 6:
- execute each swap which improves results
      immediately, and hence typically multiple swaps per iteration;
      this swapping algorithm runs in \(O(n^2)\) rather than
      \(O(n(n-k)k)\) time which is much faster for all but small \(k\). 
  
In addition, ‘FasterPAM’ uses random initialization of the
  medoids (instead of the ‘build’ phase) to avoid the
  \(O(n^2 k)\) initialization cost of the build algorithm.  In particular
  for large k, this yields a much faster algorithm, while preserving a
  similar result quality.
One may decide to use repeated random initialization by setting
  nstart > 1.