มีคนถามมาว่าช่วยทำตัวอย่างให้ดูหน่อยในการใช้ OpenCL ใน Mathematica ซึ่งผมก็บอกไปแล้วว่ามันตัวอย่างเยอะเลยในhelp Mathematica เอง 🙂 แต่ก็ยังอยากให้ผมทำให้ดู ตัวอย่างที่ทำให้ดูก็เป็นการประมาณค่า Pi โดยอาศัยที่เรียกว่า Monte carlo โดยประมาณจากการสุ่มจุดลงในสี่เหลี่ยมจัตุรัสที่มีวงกลมรัศมีครึ่งหนึ่งของความยาวด้านสี่เหลี่ยมจัตุรัส แล้วดูว่ามีกี่จุดที่อยู่ในวงกลม ซึ่งค่า Pi จะสามรถประมาณได้จาก
Pi = 4 * จำนวนจุดในวงกลม / จำนวนจุดที่สุ่มมาทั้งหมด
อันนี้เป็น code ที่ผมเขียนโดยใช้ OpenCL ใน Mathematica เพื่อหาดูว่าจุดไหนบ้างที่อยู่ในวงกลมซึ่งมีจุดศูนย์กลาง (0.5,0.5) ในสี่เหลี่ยนจัตุรัส
<< OpenCLLink`
scr = "
__kernel void myPie( __global double *X, __global double *Y,
__global mint *inC, mint np) {
int indx = get_global_id(0);
double d; //distance
if (indx < np) {
d = sqrt((X[indx]-0.5)*(X[indx]-0.5)+(Y[indx]-0.5)*(Y[indx]-0.5));
if(d <= 0.5)
inC[indx] = 1;
}
}
";
fn = OpenCLFunctionLoad[scr,
"myPie", {{_Real, _, "Input"}, {_Real, _, "Input"}, {_Integer, _,
"Output"}, _Integer}, 16 ];
np = 10^6;
x = OpenCLMemoryLoad[RandomReal[1, np]];
y = OpenCLMemoryLoad[RandomReal[1, np]];
inc = OpenCLMemoryLoad[ConstantArray[0, np]];
Print["my Pi : ", (OpenCLMemoryGet@(fn[x, y, inc, np] // First) //
Total)*4.0/np]
OpenCLMemoryUnload[x, y, inc]ทดลองรันดูสัก 100 ครั้ง โดยที่แต่ล่ะครั้งสุ่มมา 10^7 จุด แล้วดูว่ามันให้ค่าประมาณแตกต่างเท่าไหร่เมื่อเทียบกับค่าจริง (จาก Mathematica) จากรูปด้านล่าง เส้นประสีดำคือที่ได้จากการประมาณ สีแดงคือเส้นค่า Pi จาก Mathematica
ส่วนอีกอันก็เป็นการหมุนภาพ โดยที่มีจุดหมุนอยู่ที่กลางภาพครับ
scr = "
__kernel void test1( __global float * mat1, __global float * mat2,
float theta, mint width, mint height) {
int xIndex = get_global_id(0);
int yIndex = get_global_id(1);
int indx = xIndex + yIndex*width;
if (xIndex >= width || yIndex >= height)
return ;
//image centre
float x0 = width/2;
float y0 = height/2;
//relative position
int xprime = xIndex-x0;
int yprime = yIndex-y0;
float sinTheta = sin(theta);
float cosTheta = cos(theta);
int nx =(xprime*cosTheta - yprime*sinTheta + x0);
int ny =(xprime*sinTheta + yprime*cosTheta + y0);
mat2[indx] = mat1[nx + width*ny];
}
";
Clear[test];
test = OpenCLFunctionLoad[scr,
"test1", {{"Float", _, "Input"}, {"Float", _, "Output"},
"Float", _Integer, _Integer}, {16, 16},
"ShellOutputFunction" -> Print];
