作者anjackie (AN)
看板MATLAB
标题[讨论] 关於mex file的问题...
时间Tue Dec 20 13:20:42 2011
由於我的程式需要跑大量for loop（六层），
所以我用mex改到c下编写，
但是c我只有一些基础（我是物化背景），
加上Matlab里的mex file又有特别的编写要求，
在请教很多人以及查了很多资料後有了一些进展，
但是目前碰到的bug很奇怪，
我在for loop里面做矩阵运算，主要是矩阵的旋转运算，
随着我for loop给的角度不同，
而产生随角度改变的旋转矩阵，
然而我发现我给的input角度一样，
output出来的旋转array不一样，
像是这样:
for (i1 = 0; i1 < ntheta; i1++)
{
for (i2 = 0; i2 < nphi; i2++)
{
rotation(R7, phi[i2], theta[i1], 0);
printf("R7 = \n %f %f \n %f %f %f \n %f %f %f \n %f %f %f \n \n"...与下面同行
,phi[i2], theta[i1], R7[0],R7[1],R7[2],R7[3],R7[4],R7[5],R7[6],R7[7],R7[8]);
rotation(R8, phi[i2], theta[i1], 0);
printf("R8 = \n %f %f \n %f %f %f \n %f %f %f \n %f %f %f \n \n"...与下面同行
,phi[i2], theta[i1], R8[0],R8[1],R8[2],R8[3],R8[4],R8[5],R8[6],R8[7],R8[8]);
}
}
print出来的结果如下:
R7 =
0.000000 3.141593 <---- 这是给的input角度
-1.000000 0.000000 -0.000000
0.000000 1.000000 0.000000 <------这是旋转矩阵
0.000000 0.000000 -1.000000
R8 =
0.000000 3.141593 <-----同样input
-0.540302 0.841471 -0.000000
0.841471 0.540302 0.000000 <-----output旋转矩阵不同
0.000000 -0.000000 -1.000000
R7 =
15.707963 3.141593
0.857553 -0.514395 0.000000
-0.514395 -0.857553 0.000000
0.000000 -0.000000 -1.000000
R8 =
15.707963 3.141593
0.654290 0.756244 0.000000
0.756244 -0.654290 0.000000
0.000000 0.000000 -1.000000
R7 =
0.000000 0.000000
0.793480 0.608596 0.000000
-0.608596 0.793480 0.000000
0.000000 0.000000 1.000000
R8 =
0.000000 0.000000
0.701369 0.712799 0.000000
-0.712799 0.701369 0.000000
0.000000 0.000000 1.000000
R7 =
15.707963 0.000000
-0.763960 -0.645264 0.000000
0.645264 -0.763960 0.000000
0.000000 0.000000 1.000000
R8 =
15.707963 0.000000
-0.722102 0.691786 0.000000
-0.691786 -0.722102 0.000000
0.000000 0.000000 1.000000
我的旋转矩阵函式如下:
#include <stdio.h>
#include <stdio.h>
#include "mex.h"
#include <math.h>
void rotation(double *output, double a, double b, double c);
void rotation(double *output, double a, double b, double c)
{
double xx[9], yy[9], zz[9];
int i, j, k;
for (i = 0; i<9; i++)
{
xx[i] = 0;
yy[i] = 0;
zz[i] = 0;
}
xx[0] = cos(a); xx[1] = sin(a); xx[2] = 0; xx[3] = -sin(a); xx[4] = cos(a); xx[5] = 0; xx[6] = 0; xx[7] = 0; xx[8] = 1;
yy[0] = cos(b); yy[1] = 0; yy[2] = -sin(b); yy[3] = 0; yy[4] = 1; yy[5] = 0; yy[6] = sin(b); yy[7] = 0; yy[8] = cos(b);
zz[0] = cos(c); zz[1] = sin(c); zz[2] = 0; zz[3] = -sin(c); zz[4] = cos(c); zz[5] = 0; zz[6] = 0; zz[7] = 0; zz[8] = 1;
double xy[9], xyz[9];
int size = 3;
for (i=0; i<size; i++)
{
for (j = 0; j<size; j++)
{
xy[size*i+j] = 0;
xyz[size*i+j] = 0;
output[size*i+j] = 0;
}
}
for ( i = 0; i<size; i++)
{
for ( j = 0; j<size; j++)
{
for ( k = 0; k<size; k++)
{
xy[size*i+j] = xy[size*i+j] + xx[size*i+k]*yy[size*k+j];
}
}
}
for ( i = 0; i<size; i++)
{
for ( j = 0; j<size; j++)
{
for ( k = 0; k<size; k++)
{
xyz[size*i+j] = xyz[size*i+j] + xy[size*i+k]*zz[size*k+j];
output[size*i+j] = xyz[size*i+j];
}
}
}
}
落落长一大篇，因为这个程式对我很重要，
如有不足，请再推文告知，谢谢
--
※ 发信站: 批踢踢实业坊(ptt.cc)
◆ From: 140.114.44.41
※ anjackie:转录至看板 C_and_CPP 12/20 17:51
1F：→ Semisphere:你在c的测试程式码也麻烦贴上来，照这样看多半你没传 12/21 09:34
2F：→ Semisphere:址，呼叫rotation时用&R7 12/21 09:36
我原始码蛮大的，这只是一开始的旋转，想要看的测试码是指全部吗？
呼叫rotation时用&R7是指把
rotation(R7, phi[i2], theta[i1], 0);
改成rotation(&R7, phi[i2], theta[i1], 0);吗？
※ 编辑: anjackie 来自: 140.114.44.41 (12/21 11:23)
3F：→ Semisphere:ya..你rotation用mex，但却不是用matlab呼叫而是在c 12/21 12:40
4F：→ Semisphere:所以我只能猜你是在c里面测试rotation，故这测试码能贴 12/21 12:41
5F：→ Semisphere:上最好，至少也把宣告处贴出来，否则若你只是宣告错误 12/21 12:41
6F：→ Semisphere:没程式码根本没法帮你 12/21 12:42