Search code examples
cmatlabmexcomplex-numbersmatlab-compiler

Call MATLAB complex matrix function from C


I have created dll of matlab function makeSbus using matlab compiler sdk to use in C. The function takes 1 integer and 2 double matrix as input and give a complex double matrix as output. It creates sparse matrices during processing in MATLAB. This is the MATLAB source code

function Sbus = makeSbus(baseMVA, bus, gen)

on = find(gen(:, 8) > 0); %% which generators are on?
gbus = gen(on, 1);               %% what buses are they at?
    


nb = size(bus, 1);
ngon = size(on, 1);
Cg = sparse(gbus, (1:ngon)', ones(ngon, 1), nb, ngon);  


%% connection matrix.element i, j is 1 if gen on(j) at bus i is ON
    Sbus =  ( Cg * (gen(on,2 ) + 1j * gen(on, 3))   - (bus(:, 3) + 
             1j * bus(:, 4)) ) / baseMVA;

The output ofsbus is a complex double matrix of 9*1 dimension.

Using data from Case 9 here, I get the following answer from MATLAB.

https://matpower.org/docs/ref/matpower5.0/case9.html

case1='case9.m';
A=loadcase(case1);
Sbus =makeSbus(A.baseMVA,A.bus, A.gen);

Sbus= {0.7230 + 0.27030i, 1.63+ 0.0654i, 0.85 - 0.1095i, 0.0 + 0.0i, -0.90- 0.30i, 0.0+ 0.0i, -1.0 - 0.35i, 0.0 + 0.0i, -1.25- 0.50i}

I wish to call this code from C. Used matlab compiler sdk to create the necessary library files. This is my code to call the function from C

mxArray* Bus, * Gen,*var;
mxArray* sbus = NULL;

Bus = mxCreateDoubleMatrix(9, 17, mxREAL);
Gen = mxCreateDoubleMatrix(3, 25, mxREAL);
sbus = mxCreateDoubleMatrix(9, 1, mxCOMPLEX);
var = mxCreateDoubleMatrix(1, 1, mxREAL);
int data1[]={100};
double B[9][13], G[3][21];
// reads B and G matrix from a csv file


memcpy(mxGetPr(Bus), B, 9 * 13 * sizeof(double));
memcpy(mxGetPr(Gen), G, 3 * 21 * sizeof(double));
memcpy(mxGetPr(var), data1, sizeof(int));
    

if (!makeSbusInitialize()) 
   {
        fprintf(stderr, "Could not initialize the library.\n");
        return -2;
    }
    else
    {
        /* Call the library function */

        mlfMakeSbus(1, &sbus, var, Bus, Gen);
      size_t i = 0, j = 0; /* loop index variables */
      size_t r = 0, c = 0; /* variables to store the row and column 
      length of the matrix */
      double* data; /* variable to point to the double data stored 
      within the mxArray */

      /* Get the size of the matrix */
      r = mxGetM(sbus);
      c = mxGetN(sbus);
      /* Get a pointer to the double data in mxArray */
      data = mxGetPr(sbus);
       for (i = 0; i < c; i++)
       {
          for (j = 0; j < r; j++)
          {
            printf("%4.2f\t", data[j * c + i]);
          }
          printf("\n");
       }
      printf("\n");

      makeSbusTerminate();
   }

when I compile it, it shows "Error using sparse Index into matrix must be positive" .The error is when it is trying to create Cg sparse matrix inside makesbus function.

Any idea how to debug it? Also, if I get a complex matrix as my output from matlab function, how to read it from C? Is the approach I'm taking here right? Any suggestion will be appreciated.


Solution

  • Here are the problems I see with your latest edits:

    Incorrect access of sparse data (assuming sbus is sparse): Sparse matrices only physically store the non-zero elements. E.g., if the matrix size is 10x10 but only 7 elements are non-zero, then there are only 7 elements physically stored in memory, not 100 elements. If you try to access these elements with two nested for-loops like you are doing via the code data[j * c + i] where i and j ranges are based on matrix size, you will access invalid memory and crash MATLAB. You have to have special code to access the elements of a sparse matrix that takes into account the indexes of the non-zero elements that are stored in the variable. See the mxGetIr and mxGetJc functions to get access to these indexes. E.g., this code snippet would simply print the data using 1-based indexing for output:

    mwIndex *Ir, *Jc, nrow, k=0; /* k will count how far we are into the sparse data */
    Ir = mxGetIr(sbus); /* pointer to row indexes */
    Jc = mxGetJc(sbus); /* pointer to cumulative number of elements up to this column */
    for( j=0; j<n; j++ ) { /* for each column */
        nrow = Jc[j+1] - Jc[j];  /* number of non-zero row elements stored for this column */
        while(nrow--) { /* for each stored element in this column */
            printf("sbus(%d,%d) = %4.2f\n", Ir[k]+1,j+1,data[k]); /* 1-based index output */
            k++; /* increment counter */
        }
    }
    

    In the above, note that sparse index data is stored as 0-based values, even though when you print it out at the MATLAB level the indexing displayed is 1-based. The above code prints out the elements in the 1st column, then the 2nd column, etc. I.e., the elements are stored in memory in column order, so as we march through the data in memory order it gets printed out naturally in column order.

    CAVEAT: The above code assumes the sparse matrix is real. If the sparse matrix is complex, then there is the imaginary data to consider. Methods to access the real & imaginary data vary depending on whether you are using R2018a+ interleaved complex memory model or R2017b- separate complex memory model. Which version of MATLAB are you using and which memory model are you compiling with? This will affect how to write the code for element access. For the R2018a+ memory model, it is simply two elements per spot instead of one. E.g., something like this:

        printf("sbus(%d,%d) = %4.2f + %4.2f i\n", Ir[k]+1,j+1,data[2*k],data[2*k+1]); /* 1-based index output */
    

    Mismatch of matrix sizes: The mxArrays Bus and Gen are created as 9x17 and 3x25, but your C variables B and G are declared as [9][13] and [3][21]. They don't match. So during your data copying there are elements that don't get set.

    Mismatch of element order in memory: mxArray matrix elements are ordered in memory column-wise, whereas native C matrix elements are ordered in memory row-wise. In essence, they are transposes of each other in memory. You can't just memcpy the elements because the ordering is different. You would need to read the data in differently, or transpose the matrix when you do the data copy.

    Incorrect copy of int to double: data1 is type int, but the mxArray var is type double. You can't copy an int bit pattern (data1) into a double (data area of mxArray var) via memcpy ... results will be garbage. Use an assignment instead. E.g., something like this would work:

    *mxGetPr(var) = data1[0];
    

    Or you could create var directly from data1[0] like this:

    var = mxCreateDoubleScalar(data1[0]);