+ gsl_matrix_memcpy (kmeans->centers, kmeans->updated_centers);
+
+ {
+ kmeans->n = 0;
+ int i;
+ /* Step 3 */
+ gsl_vector_long_set_all (kmeans->num_elements_groups, 0.0);
+ gsl_matrix_set_all (kmeans->updated_centers, 0.0);
+ struct ccase *c;
+ struct casereader *cs = casereader_clone (reader);
+ for (; (c = casereader_read (cs)) != NULL; i++, case_unref (c))
+ {
+ int group = -1;
+ kmeans_get_nearest_group (kmeans, c, qc, &group, NULL, NULL, NULL);
+
+ for (j = 0; j < qc->n_vars; ++j)
+ {
+ const union value *val = case_data (c, qc->vars[j]);
+ if ( var_is_value_missing (qc->vars[j], val, qc->exclude))
+ continue;
+
+ double *x = gsl_matrix_ptr (kmeans->updated_centers, group, j);
+ *x += val->f;
+ }
+
+ long *n = gsl_vector_long_ptr (kmeans->num_elements_groups, group);
+ *n += qc->wv ? case_data (c, qc->wv)->f : 1.0;
+ kmeans->n++;
+
+
+ }
+ casereader_destroy (cs);
+
+
+ /* Divide the cluster sums by the number of items in each cluster */
+ for (g = 0; g < qc->ngroups; ++g)
+ {
+ for (j = 0; j < qc->n_vars; ++j)
+ {
+ long n = gsl_vector_long_get (kmeans->num_elements_groups, g);
+ double *x = gsl_matrix_ptr (kmeans->updated_centers, g, j);
+ *x /= n ;
+ }
+ }
+
+ double d = diff_matrix (kmeans->updated_centers, kmeans->centers);
+ if (d < kmeans->convergence_criteria)
+ break;
+ }