从c调用MATLAB复矩阵函数

Question

我创建了matlab函数makeSbus的dll，使用matlab编译器sdk在c中使用。该函数以1整数和2双矩阵作为输入，并给出一个复双矩阵作为输出。它在MATLAB处理过程中生成稀疏矩阵。这是MATLAB源代码

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;

sbus的输出是一个9*1维的复双矩阵。

使用案例9中的数据，我从MATLAB中得到以下答案。

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

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

全局的={0.7230 + 0.27030i，1.63 + 0.0654,0.85 - 0.1095i;0.0 + 0.0i，-0.90 - 0.30我,0.0 + 0.0,-1.0 - 0.35i;0.0 + 0.0i，-1.25 - 0.50我}

我想从c调用这个代码使用matlab编译器sdk创建必要的库文件。这是我从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.2ft", data[j * c + i]);
}
printf("n");
}
printf("n");
makeSbusTerminate();
}

当我编译时，它显示"Error使用稀疏必须positive&quot索引矩阵;.错误是当它试图在makesbus函数内创建Cg稀疏矩阵时。

知道怎么调试吗?同样，如果我得到一个复矩阵作为我从matlab函数的输出，如何从C中读取它?我现在采取的方法对吗?如有任何建议，我将不胜感激。

Answer 1

以下是我在您的最新编辑中看到的问题:

对稀疏数据的错误访问(假设sbus是稀疏的):稀疏矩阵只在物理上存储非零元素。例如，如果矩阵大小为10x10，但只有7个元素非零，则内存中物理存储的元素只有7个，而不是100个。如果您尝试使用两个嵌套的for循环访问这些元素，就像您通过代码data[j * c + i]所做的那样，其中i和j范围是基于矩阵大小的，您将访问无效的内存并崩溃MATLAB。您必须使用特殊的代码来访问稀疏矩阵的元素，该矩阵考虑存储在变量中的非零元素的索引。请参阅mxGetIr和mxGetJc函数来访问这些索引。例如，下面的代码片段将使用基于1的索引输出数据:

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.2fn", Ir[k]+1,j+1,data[k]); /* 1-based index output */
k++; /* increment counter */
}
}

在上面，请注意稀疏索引数据存储为基于0的值，即使在MATLAB级别打印它时显示的索引是基于1的。上面的代码打印出第一列中的元素，然后是第二列，以此类推。也就是说，元素按列顺序存储在内存中，因此当我们按内存顺序遍历数据时，它自然会按列顺序打印出来。

注意:上面的代码假设稀疏矩阵是实数。如果稀疏矩阵是复杂的，那么就需要考虑虚数据。方法访问真实的&虚拟数据取决于您是否使用R2018a+交错复杂内存模型或R2017b-单独复杂内存模型。你使用的是哪个版本的MATLAB，你用哪个内存模型编译?这将影响如何编写用于元素访问的代码。对于R2018a+内存模型，每个点只有两个元素，而不是一个。例如，像这样:

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

矩阵大小不匹配:mxArrays Bus和Gen被创建为9x17和3x25，但是C变量B和G被声明为[9][13]和[3][21]。它们不匹配。因此，在数据复制过程中，有些元素不会被设置。

内存中元素顺序不匹配:mxArray矩阵元素按内存列顺序排列，而原生C矩阵元素按内存行顺序排列。本质上，它们是彼此在内存中的转置。你不能仅仅因为顺序不同就记忆元素。您需要以不同的方式读取数据，或者在进行数据复制时对矩阵进行转置。

将int复制为double错误:data1是int类型，但是mxArray变量是double类型。你不能通过memcpy将int位模式(data1)复制到双(mxArray var的数据区域)…结果将是垃圾。用赋值来代替。例如:

*mxGetPr(var) = data1[0];

或者你可以直接从data1[0]创建var，像这样:

var = mxCreateDoubleScalar(data1[0]);