#c #double #matrix-multiplication #mex
#c #двойное #матрица- умножение #mex
Вопрос:
Я пытаюсь выполнить простое умножение матрицы на вектор, и по какой-то причине я получаю неправильный знак в своих результатах в паре моих умножений. Я понятия не имею, почему это происходит, любые указатели были бы весьма признательны.
Вот весь мой код, то есть функция matrix * vector и вызывающая функция (или как бы она ни называлась) из mex — я запускаю код из Matlab через mex.
#include "mex.h"
void mxv(int m, int n, double *A, double *b, double *c) {
double sum;
int i, j;
for (i = 0; i < m; i ) {
sum = 0.0;
for (j = 0; j < n; j ) {
sum = A[i * n j] * b[j];
}
c[i] = sum;
}
}
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
double *A, *b, *c;
int i, j, Am, An;
A = mxGetPr(prhs[0]);
Am = (int)mxGetM(prhs[0]);
An = (int)mxGetN(prhs[0]);
c = malloc(Am * sizeof(double));
b = mxGetPr(prhs[1]);
mxv(Am, An, A, b, c);
for (i = 0; i < Am; i )
printf("c[%d] = %1.4fn", i, c[i]);
}
Я вызываю функцию mex со следующими входными данными:
A = [-865.6634 0 0 0;
0 -17002.6822 0 0;
0 0 -1726.2421 2539.6267;
0 0 -2539.6267 -1726.2421;]
b = [-0.00153521; -0.00011165; -0.00037659; 0.00044981]
Правильный результат должен быть:
1.3290
1.8983
1.7924
0.1799
Но я получаю
c[0] = 1.3290
c[1] = 1.8983
c[2] = -0.4923
c[3] = -1.7329
Итак, первые два верны (c [0] и c[1]), но не последние два.
Я добавил в свой код кучу инструкций print, чтобы попытаться выяснить, где возникает ошибка:
#include "mex.h"
void mxv(int m, int n, double *A, double *b, double *c) {
double sum;
int i, j;
for (i = 0; i < m; i ) {
sum = 0.0;
printf("********n");
for (j = 0; j < n; j ) {
printf("A[%d][%d] = %1.10fnb[%d] = %1.10fn", i , j, A[i n * j], j, b[j]);
printf("A[%d][%d]*b[%d] = %1.10fn", i, j, j, A[i * n j] * b[j]);
sum = A[i * n j] * b[j];
printf("sum = %1.10fn", sum);
}
c[i] = sum;
}
}
void mexFunction(int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[]) {
double *A, *b, *c;
int i, j, Am, An;
A = mxGetPr(prhs[0]);
Am = (int)mxGetM(prhs[0]);
An = (int)mxGetN(prhs[0]);
printf("size(A) = (%d,%d)n", Am, An);
for (i = 0; i < Am; i ) {
for (j = 0; j < An; j ) {
printf("A[%d][%d] = %1.4fn", i, j, A[i Am * j]);
}
}
c = malloc(Am *sizeof(double));
b = mxGetPr(prhs[1]);
for (i = 0; i < Am; i ) {
printf("b[%d] = %1.4fn", i, b[i]);
}
mxv(Am, An, A, b, c);
for (i = 0; i < Am; i )
printf("c[%d] = %1.4fn", i, c[i]);
}
Я почти уверен, что получаю правильные входные данные из Matlab:
size(A) = (4,4)
A[0][0] = -865.6634
A[0][1] = 0.0000
A[0][2] = 0.0000
A[0][3] = 0.0000
A[1][0] = 0.0000
A[1][1] = -17002.6822
A[1][2] = 0.0000
A[1][3] = 0.0000
A[2][0] = 0.0000
A[2][1] = 0.0000
A[2][2] = -1726.2421
A[2][3] = 2539.6267
A[3][0] = 0.0000
A[3][1] = 0.0000
A[3][2] = -2539.6267
A[3][3] = -1726.2421
b[0] = -0.0015
b[1] = -0.0001
b[2] = -0.0004
b[3] = 0.0004
Но когда я изучаю матрично-векторное умножение, я обнаруживаю, что результаты умножения по какой-то причине в некоторых случаях имеют неправильный знак. Это то, что он выводит для i = 2:
********
A[2][0] = 0.0000000000
b[0] = -0.0015352100
A[2][0]*b[0] = -0.0000000000
sum = 0.0000000000
A[2][1] = 0.0000000000
b[1] = -0.0001116500
A[2][1]*b[1] = -0.0000000000
sum = 0.0000000000
A[2][2] = -1726.2421000000
b[2] = -0.0003765900
A[2][2]*b[2] = 0.6500855124
sum = 0.6500855124
A[2][3] = 2539.6267000000
b[3] = 0.0004498100
A[2][3]*b[3] = -1.1423494859 <- THIS SHOULD BE 1.142349... (no minus sign)
sum = -0.4922639735
********
Нечто подобное происходит для случая i = 3.
Заранее спасибо!
Комментарии:
1. Вы используете правильные индексы? Учитывая, что
A[2][3]
должно быть,2539.6267
покаA[3][2]
есть-2539.6267
. Таким образом, результат правильный, если индексыA
были переключены.2. Да, точно! Глупая ошибка от моего имени. Спасибо, что взглянули!
Ответ №1:
В:
printf("A[%d][%d] = %1.10fnb[%d] = %1.10fn", i , j, A[i n * j], j, b[j]);
printf("A[%d][%d]*b[%d] = %1.10fn", i, j, j, A[i * n j] * b[j]);
Первое printf
использует A[i n*j]
, в то время как второе использует A[i*n j]
. Это транспонированные позиции в массиве.
Комментарии:
1. Ага! Конечно, глупая ошибка! Изменение всех индексов на A [i m * j] сделало свое дело. Спасибо за помощь!