#include <mex.h>
#include <matrix.h>         //for Matlab mx and mex fuctions
#include <math.h>

//   0           0
//stopcol=AE2EA(AE)
void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray*prhs[] )
{
    //input
    double *AE;

    //output
    double *stopcol;

    //local
    int i, j, k;
    int rows;
    int cols;
    int *colBaseAddr;

    AE = mxGetPr(prhs[0]);

    rows = mxGetM(prhs[0]);
    cols = mxGetN(prhs[0]);

    plhs[0] = mxCreateDoubleScalar(0);
    stopcol = mxGetPr(plhs[0]);
    *stopcol = 0;

    colBaseAddr = (int *)mxMalloc(sizeof(int)*cols);
    for (j=0,k=0; j<cols; j++)
    {
        colBaseAddr[j] = k;
        k += rows;
    }

    for (i=0; i<rows; i++)
    {
        if (*(AE+colBaseAddr[i]+i)==0)
        {
            for (k=i+1; k<rows; k++)
            {
                if (*(AE+colBaseAddr[i]+k)==1)
                {
                    break;
                }
            }
            if (k>=rows)
            {
                *stopcol = i+1;
                return;
            }
            for (j=0; j<cols; j++)
            {
                *(AE+colBaseAddr[j]+i) = *(AE+colBaseAddr[j]+i) + *(AE+colBaseAddr[j]+k);
                if (*(AE+colBaseAddr[j]+i)==2)
                {
                    *(AE+colBaseAddr[j]+i) = 0;
                }
            }
        }
        for (k=0; k<rows; k++)
        {
            if (k!=i && *(AE+colBaseAddr[i]+k)==1)
            {
                for (j=0; j<cols; j++)
                {
                    *(AE+colBaseAddr[j]+k) = *(AE+colBaseAddr[j]+k) + *(AE+colBaseAddr[j]+i);
                    if (*(AE+colBaseAddr[j]+k)==2)
                    {
                        *(AE+colBaseAddr[j]+k) = 0;
                    }
                }
            }
        }
    }
}