bfs_mex.c

The MatlabBGL library fills a hole in Matlab s suite of algorithms. Namely, it provides a rich set o

/*
 * ==============================================================
 * bfs_mex.c The mex interface to the matlab bgl wrapper.
 *
 * David Gleich
 * 20 April 20020
 * =============================================================
 */

#include "mex.h"

#include "matlab_bgl.h"

#include <math.h>
#include <stdlib.h>

void expand_to_double(int* src, double* dst, int len, double offset)
{
    int i;
    for (i = len-1; i>=0; i--)
    {
        dst[i] = (double)src[i] + offset;
    }
}

/*
 * The mex function runs a max-flow min-cut problem.
 */
void mexFunction(int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[])
{
    int i;
    
    int mrows, ncols;
    
    int n,nz;
    
    /* sparse matrix */
    int *ia, *ja;
    
    /* source/sink */
    int u;
    
    /* output data */
    double *d, *dt, *pred;
    int *int_pred;
    int *int_d;
    int *int_dt;
    
     
    
    if (nrhs != 2) 
    {
        mexErrMsgTxt("2 inputs required.");
    }

    /* The first input must be a sparse matrix. */
    mrows = mxGetM(prhs[0]);
    ncols = mxGetN(prhs[0]);
    if (mrows != ncols ||
        !mxIsSparse(prhs[0])) 
    {
        mexErrMsgTxt("Input must be a square sparse matrix.");
    }
    
    n = mrows;
        
    /* The second input must be a scalar. */
    if (mxGetNumberOfElements(prhs[1]) > 1 || !mxIsDouble(prhs[1]))
    {
        mexErrMsgTxt("Invalid scalar.");
    }
    
    /* Get the sparse matrix */
    
    /* recall that we've transposed the matrix */
    ja = mxGetIr(prhs[0]);
    ia = mxGetJc(prhs[0]);
    
    nz = ia[n];
    
    /* Get the scalar */
    u = (int)mxGetScalar(prhs[1]);
    u = u-1;
    
    plhs[0] = mxCreateDoubleMatrix(n,1,mxREAL);
    plhs[1] = mxCreateDoubleMatrix(n,1,mxREAL);
    plhs[2] = mxCreateDoubleMatrix(1,n,mxREAL);
    
    d = mxGetPr(plhs[0]);
    dt = mxGetPr(plhs[1]);
    pred = mxGetPr(plhs[2]);
    
    int_d = (int*)d;
    int_dt = (int*)dt;
    
    for (i=0; i < n; i++)
    {
        int_d[i]=-1;
        int_dt[i]=-1;
    }
    
    int_d[u] = 0;
    
    #ifdef _DEBUG
    mexPrintf("bfs...");
    #endif 
    breadth_first_search(n, ja, ia,
        u, (int*)d, (int*)dt, (int*)pred);
    #ifdef _DEBUG
    mexPrintf("done!\n");
    #endif 
    
    int_pred = (int*)pred;
    for (i=0; i<n;i++)
    {
        if (int_pred[i] == i) { int_pred[i] = -1; }
    }
    
    expand_to_double((int*)d, d, n, 0.0);
    expand_to_double((int*)dt, dt, n, 0.0);
    expand_to_double((int*)pred, pred, n, 1.0);
    
    #ifdef _DEBUG
    mexPrintf("return\n");
    #endif 
}