The initial parameters are n_clusters, n_features, n_cases
 that describe the shape of the problem and n_outer, n_inner
 bounding the loops.

getInt proc (lowerBound in int) ⇒ (result out int)
   get result
   assert result ≥ lowerBound
end proc

nClusters const = getInt(2)
nFeatures const = getInt(1)
nCases    const = getInt(nClusters)
nOuter    const = getInt(1)
nInner    const = getInt(1)

 Next we read the feature matrix.

feature var := {{0.0} @@ nFeatures} @@ nCases
loop for caseIndex in [0,nCases)
    loop for featureIndex in [0,nFeatures)
       get feature[caseIndex][featureIndex]
    end loop
end loop

 Make a random selection of cases to serve as cluster centres.

randomAssignment proc (centre out [][]real)
   perm var []int
   loop for i in [0,nCases)
      perm[i] := i
   end loop
   loop for c in [0,nClusters)
      n const = nCases - c                 number of elements to choose from
      i const = n - 1                       last of those elements
      j const = trunc(float(n) * random())  j in [0,n)
      (perm[j], perm[i]) := (perm[i], perm[j])
      centre[c] := feature[perm[i]]
   end loop
end proc    

 The squared Euclidean distance between two feature vectors.

distanceSquared proc (x in []real, y in []real) ⇒ (d out real)
   d := 0.0
   loop for i in [0,size x)
      d +:= (x[i] - y[i])**2
   end loop
end proc

repartition proc (
   oldCentre in  [][]real,
   newCentre out [][]real,
   group      out []int       case number -> group number
) ⇒ (
   sumOfSquaredDistances out real
)
   sumOfSquaredDistances := 0.0
   tally var := {0.0} @@ nClusters
   group := {}
   newCentre := {{0.0} @@ nFeatures} @@ nClusters
   loop for caseIndex in [0,nCases)
      thisCase const = feature[caseIndex]
      bestDistanceSquared var := 0.0
      bestCluster          var := -1
      loop for c in [0,nClusters)
         d const = distanceSquared(thisCase, oldCentre[c])
         if bestCluster < 0 or d < bestDistanceSquared then
            (bestCluster, bestDistanceSquared) := (c, d)
         end if
      end if
      assert bestCluster ≥ 0
      group[caseIndex,) := {bestCluster}
      tally[bestCluster] +:= 1.0
      loop for featureIndex in [0,nFeatures)
         newCentre[bestGroup, featureIndex] +:= thisCase[featureIndex]
      end loop
      loop for cluster in [0,nClusters)
         if tally[cluster] > 0.0 then
            loop for featureIndex in [0,nFeatures)
               newCentre[cluster][featureIndex] /= tally[cluster]
            end loop
         end if
      end loop
      sumOfSquaredDistances +:= bestDistanceSquared
   end loop
end proc   

outerBestSumOfSquaredDistances var real
outerBestGroup var []int

loop for outer in [0,nOuter)
   centreOne var := {{0.0} @@ nFeatures} @@ nClusters
   centreTwo var := {{0.0} @@ nFeatures} @@ nClusters
   d1         var := 0.0
   d2         var := 0.0
   group      var []int

   randomAssignment(centreOne)
   loop for inner in [0,nInner)
      d1 := repartition(centreOne, centreTwo, group)
      d2 := repartition(centreTwo, centreOne, group)
      while d2 < d1
   end loop
   if outer = 0 or d2 < outerBestSumOfSquaredDistances then
      (outerBestGroup, outerBestSumOfSquaredDistances) := (group, d2)
   end if
end loop

print outerBestGroup