Thread: efficient covariance matrix

hello,
i have a routine to calculate the 3x3 covariance matrix of a nx3 matrix, however it's not very efficient as it uses 3 nested for loops, and with large data sets it really clogs up. i was wondering if anyone could offer any advice to speed it up at all, the code is below, cheers,

Code:

double** getcovariancemat(Object* o, Object* p, int ndims)
{
	int i, j, k;

	//malloc space for cov-matrix here
	
	Vertex* omean = (Vertex*) malloc(sizeof(Vertex));
	Vertex* pmean = mean(p);

	//number of vertices in the data object (the subsampled object)
	int pverts = p->nvertices;
	Vertex* ann[pverts];
	void* kdtree = buildkdtree(o);
	
	//create an array of pverts size filled with nearest neighbors of 
	//vertices of p contained in o
	for (i = 0; i < pverts; ++i) {
		ann[i] = nearestneighbour(&(p->vertices[i]), (struct kdtree*)kdtree);
		omean->x += (ann[i]->x);
		omean->y += (ann[i]->y);
		omean->z += (ann[i]->z);
		//printf("x: %g, y: %g, z: %g, \n",ann[i]->x, ann[i]->y, ann[i]->z);
	}

	omean->x /= pverts;
	omean->y /= pverts;
	omean->z /= pverts;
	
	//arrays for calculating the variances from the data sets
	double modelvar[pverts][3];
	double datavar[pverts][3];
	
	for (i = 0; i < pverts; ++i) {
		modelvar[i][0] = (p->vertices[i].x - pmean->x);
		modelvar[i][1] = (p->vertices[i].y - pmean->y);
		modelvar[i][2] = (p->vertices[i].z - pmean->z);
		datavar[i][0] = (ann[i]->x - omean->x);
		datavar[i][1] = (ann[i]->y - omean->y);
		datavar[i][2] = (ann[i]->z - omean->z);
	}
	
	//calculating the covariance matrix - not very efficient O(n^3)
	for (i=0; i< ndims; ++i) {
		for (j = 0; j < ndims; ++j) {
			for (k=0; k < pverts; ++k) {
				double val = (modelvar[k][i]*datavar[k][j]);
				m[i][j] += val;
			}
			m[i][j] /= (pverts);
		}
	}
	
        for (i = 0; i < pverts; ++i) {
		free(ann[i]);
	}
	free(pmean);
	free(omean);
	return m;
}