Computing tangent space basis vectors for an arbitrary mesh This is more like a share and a request than a question. I converted Eric Lengyel's code, which calculates tangents of a mesh for the purpose of texturing and normal mapping, to support SIMD. For this I used DirectXMath library from windows SDK. I don't know if this is actually faster than regular floating point version (didn't test it yet). So, my request is: please let me know if you have any suggestions as to optimize the code even further.

void CalculateTangentsAndNormals( vector&lt;Vertex&gt;&amp; verts, vector &lt;U32&gt;&amp; idx )
{
    // Computing Tangent Space Basis Vectors for an Arbitrary Mesh 
    //  http://www.terathon.com/code/tangent.html
    //

    const XMVECTOR NullVector = XMLoadFloat4A(&amp;XMFLOAT4A(0.f, 0.f, 0.f, 0.f));
    const XMVECTOR W_Null = XMLoadFloat4A(&amp;XMFLOAT4A(1.f, 1.f, 1.f, 0.f));
    const XMVECTOR signature = XMLoadFloat4A(&amp;XMFLOAT4A(1.f,-1.f,-1.f,1.f));
    XMMATRIX st(NullVector, NullVector, NullVector, NullVector);

    vector &lt;XMVECTOR&gt; vTangents, vBitangents;
    const U32 NumberOfVertices = (U32)verts.size();

    vTangents.resize( NumberOfVertices , NullVector);
    vBitangents.resize( NumberOfVertices, NullVector);

    const U32 NumberOfIndices = (U32)idx.size();
    for ( U32 i = 0; i &lt; NumberOfIndices; ++i )
    {
        const U32 i0 = idx[i];
        const U32 i1 = idx[++i];
        const U32 i2 = idx[++i];

        const XMVECTOR v0 = XMLoadFloat3A(&amp;verts[i0].pos);
        const XMVECTOR v1 = XMLoadFloat3A(&amp;verts[i1].pos);
        const XMVECTOR v2 = XMLoadFloat3A(&amp;verts[i2].pos);

        const XMVECTOR e0 = XMVectorMultiply( v1 - v0, W_Null);
        const XMVECTOR e1 = XMVectorMultiply( v2 - v0, W_Null);

        const XMMATRIX e01( e0, e1, NullVector, NullVector);
        /*      | e0.x e0.y e0.z 0 |
                | e1.x e1.y e1.z 0 |
        e01 =   | 0    0    0    0 |
                | 0    0    0    0 |
        */
        const XMVECTOR t0 = XMLoadFloat2A(&amp;verts[i0].tex);
        const XMVECTOR t1 = XMLoadFloat2A(&amp;verts[i1].tex);
        const XMVECTOR t2 = XMLoadFloat2A(&amp;verts[i2].tex);

        XMVECTOR s = XMVectorMergeXY(t1 - t0, t2 - t0); // s = (t1x - t0x, t2x - t0x, t1y - t0y, t2y - t0y)

        XMFLOAT4A d;
        XMStoreFloat4A(&amp;d, s);

        const float DetInv = 1.0f/( d.x * d.w - d.z * d.y );

        s *= DetInv;
        s = XMVectorMultiply(s, signature); // s = (sx, -sy, -sz, sw)

        const U32 *const pSrc = reinterpret_cast&lt;const U32 *const&gt;(&amp;s);
        U32* pDst = reinterpret_cast&lt;U32*&gt;(&amp;st.r[0]);
        pDst[0] = pSrc[3];  // _00 = sw
        pDst[1] = pSrc[2];  // _01 = -sz
        pDst = reinterpret_cast&lt;U32*&gt;(&amp;st.r[1]);
        pDst[0] = pSrc[1];  // _10 = -sy
        pDst[1] = pSrc[0];  // _11 = sz

        /*      | sw  -sz 0 0 |
                | -sy sx  0 0 |
        st =    | 0    0  0 0 |*DetInv
                | 0    0  0 0 |
        */

        const XMMATRIX uv( XMMatrixMultiply(st, e01) );

        vTangents[i0] += uv.r[0];
        vTangents[i1] += uv.r[0];
        vTangents[i2] += uv.r[0];

        vBitangents[i0] += uv.r[1];
        vBitangents[i1] += uv.r[1];
        vBitangents[i2] += uv.r[1];

    }
    // orthogonalize normals and tangents
    for ( U32 i = 0; i &lt; NumberOfVertices; ++i )
    {
        const XMVECTOR normal = XMVector3Normalize(XMLoadFloat3A(&amp;verts[i].normal));

        // normalize tangents and orthogonalize TBN
        XMVECTOR n0 = XMVector3Normalize( vTangents[i] - XMVector3Dot( normal, vTangents[i] ) * normal );
        //calculate handedness
        const XMVECTOR n1 = XMVector3Cross( normal, vTangents[i] );
        const float w = (XMVector3Less(XMVector3Dot(n1, vBitangents[i] ), NullVector)) ? -1.f : 1.f;

        n0 = XMVectorSetW(n0, w);
        XMStoreFloat4A(&amp;verts[i].tangent, n0);
    }
}
Here's my pass at it; without the rest of your source I can't verify it's integrity so I don't know if this will compile correctly with the rest of your code, so hopefully it can help:

// outside function since they're consts, decreases run time as no more construction each function call
const XMVECTOR NullVector = XMLoadFloat4A(&amp;XMFLOAT4A(0.f, 0.f, 0.f, 0.f));
const XMVECTOR W_Null = XMLoadFloat4A(&amp;XMFLOAT4A(1.f, 1.f, 1.f, 0.f));
const XMVECTOR signature = XMLoadFloat4A(&amp;XMFLOAT4A(1.f,-1.f,-1.f,1.f));

void CalculateTangentsAndNormals( vector&lt;Vertex&gt;&amp; verts, vector &lt;U32&gt;&amp; idx )
{
    // Computing Tangent Space Basis Vectors for an Arbitrary Mesh 
    //  http://www.terathon.com/code/tangent.html

    const U32 NumberOfVertices = (U32)verts.size();
    const U32 NumberOfIndices = (U32)idx.size();
    U32 i = 0;

    XMMATRIX st(NullVector, NullVector, NullVector, NullVector);

    // explicit construction
    vector &lt;XMVECTOR&gt; vTangents( NumberOfVertices , NullVector);
    vector &lt;XMVECTOR&gt; vBitangents( NumberOfVertices, NullVector);

    {
        U32 i0, i1, i2;
        XMVECTOR v0, v1, v2;
        XMVECTOR t0;
        XMFLOAT4A d;
        const U32* pSrc = 0;
        U32* pDst = 0;
        while (i &lt; NumberOfIndices) {
            i0 = idx[i];
            i1 = idx[++i];
            i2 = idx[++i];

            v0 = XMLoadFloat3A(&amp;verts[i0].pos);
            v1 = XMLoadFloat3A(&amp;verts[i1].pos);
            v2 = XMLoadFloat3A(&amp;verts[i2].pos);

            // e0, e1, NullVector, NullVector
            const XMMATRIX e01( XMVectorMultiply( v1 - v0, W_Null),
                                XMVectorMultiply( v2 - v0, W_Null),
                                NullVector,
                                NullVector);
            /*      | e0.x e0.y e0.z 0 |
                    | e1.x e1.y e1.z 0 |
            e01 =   | 0    0    0    0 |
                    | 0    0    0    0 |
            */
            t0 = XMLoadFloat2A(&amp;verts[i0].tex);

            // s = (t1x - t0x, t2x - t0x, t1y - t0y, t2y - t0y)
            XMVECTOR s = XMVectorMergeXY((XMLoadFloat2A(&amp;verts[i1].tex) - t0),
                                         (XMLoadFloat2A(&amp;verts[i2].tex) - t0));
            XMStoreFloat4A(&amp;d, s);
            s *= (1.0f/( d.x * d.w - d.z * d.y )); //DetInv;
            s = XMVectorMultiply(s, signature); // s = (sx, -sy, -sz, sw)

            pSrc = reinterpret_cast&lt;const U32*&gt;(&amp;s);
            pDst = reinterpret_cast&lt;U32*&gt;(&amp;st.r[0]);
            pDst[0] = pSrc[3];  // _00 = sw
            pDst[1] = pSrc[2];  // _01 = -sz
            pDst = reinterpret_cast&lt;U32*&gt;(&amp;st.r[1]);
            pDst[0] = pSrc[1];  // _10 = -sy
            pDst[1] = pSrc[0];  // _11 = sz

            /*      | sw  -sz 0 0 |
                    | -sy sx  0 0 |
            st =    | 0    0  0 0 |*DetInv
                    | 0    0  0 0 |
            */

            const XMMATRIX uv( XMMatrixMultiply(st, e01) );

            vTangents[i0] += uv.r[0];
            vTangents[i1] += uv.r[0];
            vTangents[i2] += uv.r[0];

            vBitangents[i0] += uv.r[1];
            vBitangents[i1] += uv.r[1];
            vBitangents[i2] += uv.r[1];

            ++i;
        }
    } // end scope for temp vars
    i = 0;
    XMVECTOR normal, n0, n1;
    float w;
    while (i &lt; NumberOfVertices) {
        normal = XMVector3Normalize(XMLoadFloat3A(&amp;verts[i].normal));
        // normalize tangents and orthogonalize TBN
        n0 = XMVector3Normalize( vTangents[i] - XMVector3Dot( normal, vTangents[i] ) * normal );
        //calculate handedness
        n1 = XMVector3Cross( normal, vTangents[i] );
        w = (XMVector3Less(XMVector3Dot(n1, vBitangents[i] ), NullVector)) ? -1.f : 1.f;
        n0 = XMVectorSetW(n0, w);
        XMStoreFloat4A(&amp;verts[i].tangent, n0);
        ++i;
    }
}