[討論] 關於mex file的問題...

[anjackie (AN)]

由於我的程式需要跑大量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
→ Semisphere:你在c的測試程式碼也麻煩貼上來，照這樣看多半你沒傳
→ Semisphere:址，呼叫rotation時用&R7
我原始碼蠻大的，這只是一開始的旋轉，想要看的測試碼是指全部嗎？
呼叫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)
→ Semisphere:ya..你rotation用mex，但卻不是用matlab呼叫而是在c
→ Semisphere:所以我只能猜你是在c裡面測試rotation，故這測試碼能貼
→ Semisphere:上最好，至少也把宣告處貼出來，否則若你只是宣告錯誤
→ Semisphere:沒程式碼根本沒法幫你