The Gustafson–Kessel’s c-means algorithm

Input: the set 𝒜 ⊂ (ℝ^n), initial centers (Subscript[z, 1]),⃛,(Subscript[z, k]) ∈ (ℝ^n);
Output: centers (Subscript[c, 1]),⃛,(Subscript[c, k]) ∈ (ℝ^n), membership matrix U, number of iterations;
Submodules: mu[], Fz[]

cmeans[A_,z_,Ind_:0]:=Module[{q=2,k=Length[z],m=Length[A],F1,F2, cen,eps=.005,It=20},
					(* Initial approximation *)
kk=1;
cen=z//N;
U1=Table[mu[A,cen,i,j],{i,m},{j,k}];  
F1=Fz[A,cen,U1];
If[Ind≠0,slcen=ListPlot[cen,PlotStyle->{Red,PointSize[.04],Opacity[.7]}];
sl2=Show[slA,slcen];
Print[GraphicsGrid[{{slA,sl2}},ImageSize->Medium]]];
						(* While loop *)
While[
cenN=Table[Sum[mu[A,cen,i,j]^q A[[i]],{i,m}]/Sum[mu[A,cen,i,j]^q,{i,m}],{j,k}];
If[Ind≠0,slcen=ListPlot[cenN,PlotMarkers->{"✶",20},PlotStyle->{Orange,PointSize[.04],Opacity[.7]}];
sl1=sl2; sl2=Show[slA,slcen];
Print[GraphicsGrid[{{sl1,sl2}},ImageSize->Medium]]];
U2=Table[mu[A,cenN,i,j],{i,m},{j,k}]; 
 DelU=U2-U1; 
U1=U2;
F2=Fz[A,cenN,U1];
If[Ind≠0,Print["It-",kk,": centri=",cenN,"; Fz = ",F2,"; accuracy = ",Norm[DelU]  ]];
Norm[DelU]>eps (*Abs[(F1-F2)/F2]>eps*)&& kk<It,F1=F2;
kk=kk+1;cen=cenN
 ];
{cen,U2,kk,(F1-F2)/F2}
]