#include "mex.h"

#include "matrix.h"

#include "papi.h"

static long long accum_error = 0;

static long long start_time = 0;

void mexFunction(int nlhs, mxArray *plhs[],

                 int nrhs, const mxArray *prhs[]) {

  float real_time, proc_time, rate;

  double *x;

  unsigned int mrows, ncols;

  int result;

  unsigned int flop_events[2];

  long long ins = 0, flop_values[2];

  long long elapsed_time;

  /* Check for proper number of arguments. */

    if(nrhs > 1) {

      mexErrMsgTxt("This function expects one optional input.");

    } else if(nlhs > 2) {

	mexErrMsgTxt("This function produces 1 or 2 outputs: [ops, mflops].");

    }

    /* The input must be a noncomplex scalar double.*/

    if(nrhs == 1) {

      mrows = mxGetM(prhs[0]);

      ncols = mxGetN(prhs[0]);

      if(!mxIsDouble(prhs[0]) || mxIsComplex(prhs[0]) || !(mrows == 1 && ncols == 1)) {

        mexErrMsgTxt("Input must be a noncomplex scalar double.");

      }

      /* Assign a pointer to the input. */

      x = mxGetPr(prhs[0]);

      /* if input is 0, reset the counters by calling PAPI_stop_counters with 0 values */

      if(*x == 0) {

         if (start_time == 0) {

            PAPI_stop_counters(NULL, 0);

            accum_error = 0;

         } else {

            start_time = 0;

            PAPI_stop_counters(flop_values, 2);

         }

      }

    }

   if(result = PAPI_event_name_to_code("EMON_SSE_SSE2_COMP_INST_RETIRED_PACKED_DOUBLE", &(flop_events[0])) < PAPI_OK) {

      if(result = PAPI_flops( &real_time, &proc_time, &ins, &rate)

         mexPrintf("Error code: %d\n", result);

	      mexErrMsgTxt("Error getting flops.");

      }

   } else {

      if(start_time == 0) {

         flop_events[1] = PAPI_FP_OPS;

         start_time = PAPI_get_real_usec();

         if((result = PAPI_start_counters(flop_events, 2)) < PAPI_OK) {

            mexPrintf("Error code: %d\n", result);

            mexErrMsgTxt("Error getting flops.");

         } else {

            ins = 0;

            rate = 0;

         }

      } else {

         if((result = PAPI_read_counters(flop_values, 2)) < PAPI_OK) {

            mexPrintf("%d\n", result);

            mexErrMsgTxt("Error reading the running counters.");

         } else {

            elapsed_time = PAPI_get_real_usec() - start_time;

            ins = (2*flop_values[0])+flop_values[1];

            rate = ((float)ins)/((float)elapsed_time);

         }

      }

   }

/*    mexPrintf("real: %f, proc: %f, rate: %f, ins: %lld\n", real_time, proc_time, rate, ins); */

    if(nlhs > 0) {

      plhs[0] = mxCreateScalarDouble((double)(ins - accum_error));

      /* this call adds 7 fp instructions to the total */

      /* but apparently not on Pentium M with Matlab 7.0.4 */

//      accum_error += 7;

      if(nlhs == 2) {

        plhs[1] = mxCreateScalarDouble((double)rate);

      /* the second call adds 4 fp instructions to the total */

      /* but apparently not on Pentium M with Matlab 7.0.4 */

//        accum_error += 4;

      }

    }

}