/* * * mexCCACollectdata.c *      prepare data for cca.m to form SDP * * by feisha@cis.upenn.edu */  #include "mex.h" #include "matrix.h"#include <stdlib.h>#include <float.h>#include <string.h> #include <math.h> /* the computation engine */void collectdata(double *x, double *y, int* edgerow, int *edgecol, int *relative, int *nv, int *vidx, int D, int d, int n, int ks, double *a, double *b, double *g);/* auxiliary functions */void sanity_check(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]);/* the gateway */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){    int n, ks=1, D, d;    double *x, *y, *a, *b, *g;    int *neighbors;        /* sanity_check(nlhs, plhs, nrhs, prhs); */        n = mxGetN(prhs[0]);    D = mxGetM(prhs[0]);    d = mxGetM(prhs[1]);    /*printf("%d data points, reducing from dimension %d ---> dimension %d, using \%d nearest neighbors\n", n, D, d, ks); */        plhs[0] = mxCreateDoubleMatrix(d*d, d*d, mxREAL);    plhs[1] = mxCreateDoubleMatrix(d*d, n, mxREAL);    plhs[2] = mxCreateDoubleMatrix(n,1, mxREAL);        collectdata( (double*)mxGetData(prhs[0]),                 (double*)mxGetData(prhs[1]),                 (int*)mxGetData(prhs[2]),                 (int*)mxGetData(prhs[3]),                 (int*)mxGetData(prhs[4]),                 (int*)mxGetData(prhs[5]),                 (int*)mxGetData(prhs[6]),                 D, d, n, ks,                  (double*)mxGetData(plhs[0]),                 (double*)mxGetData(plhs[1]),                 (double*)mxGetData(plhs[2]));    return;}/* check input arguments */void sanity_check(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {    int n1, n2, n3;    /* Check for proper number of input and output arguments. */        if (nrhs != 4) {        mexErrMsgTxt("Four input argument required.");    }     if (nlhs != 3) {        mexErrMsgTxt("Three output arguments required.");    }    /* Check data type of input argument. */    /*if (!(mxIsDouble(prhs[0])) || !(mxIsDouble(prhs[1])) || !(mxIsDouble(prhs[2]))) {        mexErrMsgTxt("Input array must be of type double.");    }*/        /* Check dimensions */    if(mxGetNumberOfDimensions(prhs[0])!=2 || mxGetNumberOfDimensions(prhs[1])!=2 || mxGetNumberOfDimensions(prhs[2])!=2) {        mexErrMsgTxt("Input arrays must be two-dimension arrays");    }    /* matching dimensions */    n1 = mxGetN(prhs[0]);    n2 = mxGetN(prhs[1]);    n3 = mxGetN(prhs[2]);    if (n1 != n2 || n1 !=n3) {        mexErrMsgTxt("Input arrays must have same number of columns!");    }    return;}/* computing */double vecdot(double* v1, int n) {  /* v1'*v2 */    double x = 0;    int i;    for(i = 0; i < n; i++) {        x += v1[i]*v1[i];    }    return x;}void vecsub(double* v1,double* v2, double* v3, int n) { /* v3=v1-v2 */    int incx=1;    double none = -1;    dcopy_(&n, v1, &incx,v3, &incx);    daxpy_(&n, &none, v2, &incx, v3, &incx);    }void vecout(double* v1, double* v3, int d) { /* v3 = v1*v2' */    double one = 1.0;    int inca = 1;    memset(v3, 0, sizeof(double)*d*d);    dger_(&d , &d, &one, v1, &inca,v1, &inca, v3, &d); }void makevecsym(double *v, int d) { /* assume V is a dxd matrix, return (V+V')/2 */    int row, col, i;    double temp;    for(col=0; col < d; col++) {        for(row=col; row < d; row++) {            temp = (v[row*d+col] + v[col*d+row])/2;            v[row*d+col] = v[col*d+row] = temp;        }    }}void updateA(double *a, double *v, double alpha, int d) {  /* update A */    int inca = 1;    char *uplo= "U";    dspr_(uplo, &d , &alpha, v, &inca, a);}void updateB(double * b, double alpha, double *v, int d) { /* update B(:,i) */    int i;    int incx = 1;    daxpy_(&d, &alpha, v, &incx, b, &incx);}/* transform upper triangular matrix stored in A as full matrix */void recoverA(int D, double *A, double *tempA){  int l = 0, idx, jdx,s, incx=1, incy=1;  int n = D*(D+1)/2;  for(jdx=0; jdx <D; jdx++) {     for(idx=0; idx <=jdx; idx++) {    	A[jdx*D+idx] = tempA[l];        l++;          }  }  for(jdx=0; jdx <D; jdx++)     for(idx=jdx; idx <D; idx++) {        A[jdx*D+idx] = A[idx*D+jdx];    } }void collectdata(double *x, double *y, int* edgesrow, int *edgescol, int *relative, int *nv, int *vidx, int D, int d, int n, int ks, double *a, double *b, double *g) {    double *diffx, *diffy, *yij, *nid;    double gij;    int i, j, k, iEdge, iVertex,itemp, ii;    int nn_start, nn_end, nn;    double *diffxjk, *diffyjk;    double gjk, *yjk, *tempA;    int relflg = *relative;    /* clear data area */    memset(a, 0, sizeof(double)*d*d*d*d);    memset(b, 0, sizeof(double)*d*d*n);    memset(g, 0, sizeof(double)*n);    /* temporary working space */    diffx = (double*)calloc(D, sizeof(double));    diffy = (double*)calloc(d, sizeof(double));    yij = (double*)calloc(d*d, sizeof(double));            tempA = (double*)calloc((d*d)*(d*d+1)/2, sizeof(double));        if(!diffx || !diffy || !yij || !tempA) {        mexErrMsgTxt("Out of memory..cannot allocate working space..");        return;    }    iEdge = 0;    iVertex = 0;    if(relflg==0) {            for(i=0; i <n ; i++) {             /* figure out the nearest neighbors */            nn_start = edgescol[i]; nn_end = edgescol[i+1]-1;            if (nn_end >= nn_start) {                for(nn=nn_start; nn<=nn_end; nn++) {                    j = edgesrow[nn];                                        /* compute diffx and diffy */                    vecsub(x+i*D, x+j*D, diffx, D);                    gij = vecdot(diffx, D);                                    vecsub(y+i*d, y+j*d, diffy, d);                    vecout(diffy,  yij, d);                                /* update A, b */                    updateA(tempA, yij, 1.0*nv[iEdge],d*d);                    for(itemp=0; itemp < nv[iEdge]; itemp++) {                        ii = vidx[iVertex]-1;                        updateB(b+ii*d*d, gij, yij, d*d);                        g[ii] = g[ii]+gij*gij;                        iVertex++;                    }                                        iEdge++;                }               }        }        printf("Iedge :%d, ivertex: %d\n", iEdge, iVertex);    } else {        /* reproducing the code above to some degree so that we don't have to          do if statement inside cache-intensive loop */        for(i=0; i <n ; i++) {             /* figure out the nearest neighbors */            nn_start = edgescol[i]; nn_end = edgescol[i+1]-1;            if (nn_end >= nn_start) {                for(nn=nn_start; nn<=nn_end; nn++) {                    j = edgesrow[nn];                                /* compute diffx and diffy */                    vecsub(x+i*D, x+j*D, diffx, D);                    gij = vecdot(diffx, D);                                    vecsub(y+i*d, y+j*d, diffy, d);                    vecout(diffy,  yij, d);                                /* update A, b */                    updateA(tempA, yij, 1*nv[iEdge]/(gij*gij), d*d);                    for(itemp=0; itemp < nv[iEdge]; itemp++) {                        ii = vidx[iVertex]-1;                        updateB(b+ii*d*d, 1/gij, yij, d*d);                        g[ii] = g[ii]+1;                        iVertex++;                    }                    iEdge++;                }               }        }    }    /* make A symmetric */    recoverA(d*d, a, tempA); }