+covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c)
+{
+ size_t i, j;
+ const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
+
+ assert (cov->passes == 2);
+ assert (cov->state >= 1);
+
+ if (! cov->pass_two_first_case_seen)
+ {
+ assert (cov->state == 1);
+ cov->state = 2;
+
+ /* Divide the means by the number of samples */
+ for (i = 0; i < cov->n_vars; ++i)
+ {
+ for (j = 0; j < cov->n_vars; ++j)
+ {
+ double *x = gsl_matrix_ptr (cov->moments[MOMENT_MEAN], i, j);
+ *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
+ }
+ }
+ }
+
+ for (i = 0 ; i < cov->n_vars; ++i)
+ {
+ const union value *val1 = case_data (c, cov->vars[i]);
+
+ if ( var_is_value_missing (cov->vars[i], val1, cov->exclude))
+ continue;
+
+ for (j = 0 ; j < cov->n_vars; ++j)
+ {
+ int idx;
+ double ss ;
+ const union value *val2 = case_data (c, cov->vars[j]);
+
+ const double s = pow2 (val1->f - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j)) * weight;
+
+ if ( var_is_value_missing (cov->vars[j], val2, cov->exclude))
+ continue;
+
+ {
+ double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
+ *x += s;
+ }
+
+ ss =
+ (val1->f - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
+ *
+ (val2->f - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
+ * weight
+ ;
+
+ idx = cm_idx (cov, i, j);
+ if (idx >= 0)
+ {
+ cov->cm [idx] += ss;
+ }
+
+ }
+ }
+
+ cov->pass_two_first_case_seen = true;
+}
+
+
+/* Call this function for every case in the data set.
+ After all cases have been passed, call covariance_calculate
+ */
+void