C mex file - sort and store index

I am trying to sort an integer index based on some real number values. The vector of real numbers is extracted from a mxArray input into my mex function. And I want to loop over the process for each row of the matrix and store the final index set in another matrix that is output for the mex file.

Utility(MxN)

RealSort (1,N)

Index (1,N)

Output (MxN)

I've thought about using pointers, but I cannot figure out how to preserve the index and not the address.

I've also tried a global "base_arr = row of Utility" and pass "Index" to qsort, but that isn't coming up with the right rank of index.

/* Include files */
#include "math.h"
#include "mex.h"

#include <stdio.h>
#include <stdlib.h>

/* Checks whether inputs are matrices */
#define IS_REAL_2D_FULL_DOUBLE(P) (!mxIsComplex(P) && \
mxGetNumberOfDimensions(P) == 2 && !mxIsSparse(P) && mxIsDouble(P))
#define IS_REAL_SCALAR(P) (IS_REAL_2D_FULL_DOUBLE(P) && mxGetNumberOfElements(P) == 1)

/*************************************************/
/* SORTING */

/* Sort elements with qsort */
int comp(const void * elem1, const void * elem2) 
{
double f = *((double*)elem1);
double s = *((double*)elem2);
    if (f > s) return -1;
    if (f < s) return 1;
    return 0;
}

double *base_arr;

static int compIndex(const void *a, const void *b)
{
  int aa = *((int *) a), bb = *((int *) b);
  if (base_arr[aa] < base_arr[bb])
    return -1;
  if (base_arr[aa] == base_arr[bb])
    return 0;
  if (base_arr[aa] > base_arr[bb])
    return 1;
}

void sort(double *y, int size)
{
  int temp;
  for (int i = 0; i < size; i ++){
    for(int j = i + 1; j < size; j++){
      if(y[j] < y[i]){
        temp = y[j];
        y[j] = y[i];
        y[i] = temp;
      }
    }
  }
}

/*************************************************/

void assignVectorWid(double *y, double *z, int A, int B, int a)
{  
  for (int b = 0; b < B; b++)
  {
    z[b] = y[b];
/*    z[b] = y[a + A*b];*/
  }
}
void assignVectorLen(double *y, double *z, int A, int B, int b)
{  
  for (int a = 0; a < A; a++)
  {
    z[a] = y[a + A*b];
  }
}

/*************************************************/

/* MAIN */

/*************************************************/

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
  /* Define output matrix */

  #define Bset     plhs[0]
  #define SortUtil   plhs[1]
  #define Rank     plhs[2]

  /* Define input matrices */

  #define Utility   prhs[0]

  /* Check the number of inputs is correct (Aoffer,Baccept,QbarA) */
  if(nrhs < 1)
    mexErrMsgTxt("Too few input arguments.");
  else if(nrhs > 1)
    mexErrMsgTxt("Too many input arguments.");

  if(!IS_REAL_2D_FULL_DOUBLE(Utility))
    mexErrMsgTxt("Utility must be real 2D full double array.");

  /* Dimensions and pointers for inputs */

  int A = mxGetM(Utility);
  int B = mxGetN(Utility);
  int nb = A*B;
  double *Utility_pt = mxGetPr(Utility);

  /* Output matrix */

  Bset = mxCreateDoubleMatrix(A,B,mxREAL);
  double *Bset_pt = mxGetPr(Bset);
  SortUtil = mxCreateDoubleMatrix(1,B,mxREAL);
  double *SortUtil_pt = mxGetPr(SortUtil);
  Rank = mxCreateDoubleMatrix(1,B,mxREAL);
  double *Rank_pt = mxGetPr(Rank);

  /* Generate index to use for sorting */

  int  *idx = malloc (sizeof (int) * B);
  for (int b = 0; b < B; b++)
  {
    idx[b] = b + 1;
    printf("index %d\n", idx[b]);
  }

/***************************************************************/

/* Begin Sorting

/***************************************************************/

  int a = 0;
  base_arr = (double *)mxGetData(Utility);
  int *idxSort = idx;
  qsort (idx, B, sizeof(*idx), compIndex);
      assignVectorWid(idx,SortUtil_pt,A,B,a);

Any help would be greatly appreciated! Thank you!

Answer 1

[James Tursa]

You compIndex routine does not appear to index into the non-contiguous row data properly. You create an index array idx that has the numbers 1 ... N, all well and good. But then inside the compIndex routine you assume simple linear indexing into the base_arr using these indexes (e.g., base_arr[aa] < base_arr[bb]). That would only work if the data in base_arr that you are accessing was contiguous in memory (e.g., columns), but the data is NOT contiguous in memory because you are attempting to sort rows. Also, you are using 1-based indexes to index into base_arr, but you need to use 0-based indexes in C. So this will not work.

You need to fix your indexing in compIndex if you really want to sort the rows. E.g., notionally

base_arr[A*(aa-1)] < base_arr[A*(bb-1)]

The A factor is to skip values by according to the number of elements in each column. The -1 is to convert the 1-based index into a 0-based index. The above notional code would only work for the first row. For other rows you would need to alter the value of base_arr, or add in the row number inside the [ ].

Finally, remember that to store the result in an output array by rows, you will have the same non-contiguous indexing issue to deal with.

The whole thing would be a lot simpler if you passed in the transpose and did the sorting by columns instead of rows. Then you could pass in the column pointers directly to qsort, the linear indexing you are attempting to use would work, and there would be no need to copy the results over to an output array since qsort could sort the output array columns directly.

Is there a reason you are using malloc instead of mxMalloc?

[abv]

Amazing! Thank you so much as it is definitely working on the row and now I can expand to the matrix and perhaps adapt what you mention about the transpose to compare. Your explanation is clear too.

I had used malloc to test the difference, and I changed it to mxMalloc based on your comment. This is my first mex adventure (and really C adventure, as well), and I was trying to map what I was reading about arrays to mxArrays and see if I'd have more success.

Thank you again!