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; }