高效率的3D图形数学库(2)---SSE与矩阵相乘

这次将介绍SSE扩展指令集，以及矩阵乘法的优化，不喜汇编者请发送WM_CLOSE消息！！闲话就不说了，SSE指令的历史到处都是，主要说说我对指令集的原理、作用和用法的理解。SSE指令集的最大有点就是能够 4个float并行运算，与其说这恰好符合图形算法，到不如说就是为图形编程设计的。比如说，两个向量相加，x1+x2, y1+y2, z1+z2, w1+w2。用了4条加法指令，不如来一条省事。这就是SSE能为我们做的。这里要介绍8个寄存器：xmm0, xmm1, xmm2...xmm6, xmm7，这些寄存器都是128位的，每个寄存器都可以存放4个float，所以，x1, y1, z1, w1这4个float型可以放在xmm0中，而x2, y2, z2, w2可以放在xmm1中，于是 xmm0 += xmm1，结果就保存在xmm0寄存器了！下面是代码:struct Vector{float x, y, z, w;void Add(const Vector* pIn){_asm{mov eax, pIn; // 这里其实应该是 mov eax, dword ptr[pIn];mov ecx, this;movups xmm0, [ecx]; // 把this的xyzw放入xmm0movups xmm1, [eax]; // 把pIn的xyzw放入xmm1addps xmm0, xmm1; // xmm0 += xmm1movups [ecx], xmm0; // 把xmm0的值给thismov [ecx+12], 3F800000h // 别忘了给w = 1.0f}}}其实也就这么简单，movups是把内存中的向量放入寄存器，或者把寄存器中的向量放回内存，总之是一次移动128个位的一条指令,这个指令比普通的MOV指令要慢70%左右，比浮点乘法要慢的多了，大量的时间花在了movups上。所以说，简单的算法不值得去用SSE指令。这个向量加法指令很可能会比下面的算法更慢：void Add(const Vector* pIn){x += pIn->x;y += pIn->y;z += pIn->z;w += pIn->w;}那么，如果我要让一个向量的xyzw这4个float同时去加上同一个float值怎么办呢？比方说，我要让在xmm0中的xyzw同时加上float型变量h，就会像下面这样去做：float h = ...;_asm{movss xmm1, hshufps xmm1, xmm1, 0addps xmm0, xmm1}好了，出现2个新的指令： movss和shufps。movss是将一个32位的值移动到一个128位寄存器的低32位中。也就是说，这时xmm1中的4个32位区块中只有第0个区块是存放了h的值。这时，我想让其他3个区块也存放同样的值，以便对xmm0进行并行加法处理，于是用到了下面的指令————shufps是用来将128位寄存器中4个区块的数值进行相互调换、拷贝或覆盖，就像洗牌一样，所以叫做洗牌指令。我们现看该指令的第3个操作数是 0，该操作数是8位的，0 = 00 00 00 00，你看，我把8位数分成了4份，每一份都表示了一个寄存器区间。00就是第0个区块01就是第1个区块10就是第2个区块11就是第3个区块shufps dest, src, 00 00 00 00我们将根据第3个操作数，从src中选取区间，并把该区块中的数给拖到dest相应的区块内。这里的操作数4个都是0，所以，dest的4个区块中存放的都是src中第0个区间的值。而dest和src是同一个寄存器的时候，就会是自我洗牌。于是xmm1中的4个区块存放的都是xmm1第0个区块的值。下面的示例会帮助你更好的理解这个问题：dest = w1, z1, y1, x1src = w2, z2, y2, x2经过该指令： shufps dest, src, 01 11 00 10得到该结果： dest = y2, w2, x2, z2看下面的会更清晰：dest第3区块 _14 = pIn1->_11*pIn2._14 + pIn1->_12*pIn2._24 + pIn1->_13*pIn2._34 + pIn1->_14*pIn2._44;pOut->_21 = pIn1->_21*pIn2._11 + pIn1->_22*pIn2._21 + pIn1->_23*pIn2._31 + pIn1->_24*pIn2._41;pOut->_22 = pIn1->_21*pIn2._12 + pIn1->_22*pIn2._22 + pIn1->_23*pIn2._32 + pIn1->_24*pIn2._42;pOut->_23 = pIn1->_21*pIn2._13 + pIn1->_22*pIn2._23 + pIn1->_23*pIn2._33 + pIn1->_24*pIn2._43;pOut->_24 = pIn1->_21*pIn2._14 + pIn1->_22*pIn2._24 + pIn1->_23*pIn2._34 + pIn1->_24*pIn2._44;pOut->_31 = pIn1->_31*pIn2._11 + pIn1->_32*pIn2._21 + pIn1->_33*pIn2._31 + pIn1->_34*pIn2._41;pOut->_32 = pIn1->_31*pIn2._12 + pIn1->_32*pIn2._22 + pIn1->_33*pIn2._32 + pIn1->_34*pIn2._42;pOut->_33 = pIn1->_31*pIn2._13 + pIn1->_32*pIn2._23 + pIn1->_33*pIn2._33 + pIn1->_34*pIn2._43;pOut->_34 = pIn1->_31*pIn2._14 + pIn1->_32*pIn2._24 + pIn1->_33*pIn2._34 + pIn1->_34*pIn2._44;pOut->_41 = pIn1->_41*pIn2._11 + pIn1->_42*pIn2._21 + pIn1->_43*pIn2._31 + pIn1->_44*pIn2._41;pOut->_42 = pIn1->_41*pIn2._12 + pIn1->_42*pIn2._22 + pIn1->_43*pIn2._32 + pIn1->_44*pIn2._42;pOut->_43 = pIn1->_41*pIn2._13 + pIn1->_42*pIn2._23 + pIn1->_43*pIn2._33 + pIn1->_44*pIn2._43;pOut->_44 = pIn1->_41*pIn2._14 + pIn1->_42*pIn2._24 + pIn1->_43*pIn2._34 + pIn1->_44*pIn2._44;}else{_asm{mov edx, pIn2; // 这时保存的是pIn2movups xmm4, [edx]; //pIn2的第1行movups xmm5, [edx+16]; //pIn2的第2行movups xmm6, [edx+32]; //pIn2的第3行movups xmm7, [edx+48]; //pIn2的第4行mov eax, pIn1; // 这时保存的是pIn1mov edx, pOut;mov ecx, 4; // 循环4次LOOPIT: // 开始循环movss xmm0, [eax]; xmm0 = pIn1->xshufps xmm0, xmm0, 0; 洗牌xmm0 = pIn1->x, pIn1->x, pIn1->x, pIn1->xmulps xmm0, xmm4;movss xmm1, [eax+4]; xmm1 = pIn1->yshufps xmm1, xmm1, 0; 洗牌xmm1 = pIn1->y, pIn1->y, pIn1->y, pIn1->ymulps xmm1, xmm5;movss xmm2, [eax+8]; xmm2 = pIn1->zshufps xmm2, xmm2, 0; 洗牌xmm2 = pIn1->z, pIn1->z, pIn1->z, pIn1->zmulps xmm2, xmm6;movss xmm3, [eax+12]; xmm3 = pIn1->wshufps xmm3, xmm3, 0; 洗牌xmm3 = pIn1->w, pIn1->w, pIn1->w, pIn1->wmulps xmm3, xmm7;addps xmm0, xmm1;addps xmm2, xmm3;addps xmm0, xmm2; 最终结果行保存在xmm0movups [edx], xmm0; 将结果保存到pOut中add edx, 16;add eax, 16; 作为变址用loop LOOPIT;}}}这里的汇编对照上面的一般算法就很容易理解。下回我会说说矩阵类的关键算法。To Be Continued.....