This section of the archives stores flipcode's complete Developer Toolbox collection, featuring a variety of mini-articles and source code contributions from our readers.

 
  Matrix Mult Tests
  Submitted by



I was at an AMD performance talk this past week, and they were showing how one of their perf tools. To demonstrate, they had 3 different matrix multiplication implementations which they ran several million times. One implementation was AMD's 3dnow implementation, another was C code, and the other D3DXMatrixMultiply. Out of the 3, obviously the 3dnow version was the fastest. However, the C code was faster than what I thought would be optimized assembly in D3DXMatrixMultiply.

I was skeptical so I wrote up an inline assembly implementation and a C version and compared them vs D3DXMatrixMultiply. Running the test, I got the following results (measuing using QueryPerformanceCounter over 5 million matrix mult operations): Asm: 1923883
C : 3669129
D3DX: 3755883 Which lead me to believe that the D3DX implementation is only vanilla C... which is crap compared to assembly =P I sent my program to the work computer the next day to test them vs the latest build of D3DX, dx8.1 (i'm currently working for microsoft/directx), and the D3DX implementation was now around 15% faster vs assembly on the same machine. I traced into D3DXMatrixMultiply and the code was completely different from what i saw on my machine at home (dx8.0a); 3dnow ops were being used. I found out later that day that the D3DX team has been playing around with processor specific optimizations, and D3DX detects your proc and uses functions specific to the type detected. (which means there are SSE/SSE2 optimizations as well if u have p3/4) Anyways, this all fits into COTD cause I have the assembly implementation attached (along with the tests), in case anyone wanted it.

Currently browsing [multmatrixasm.zip] (4,499 bytes) - [matmult.cpp] - (16,204 bytes)

// matmult.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include <memory.h>
#include <iostream>
#include <windows.h>
#include <D3dx8math.h>
#include <stdio.h>
#include <assert.h>

using namespace std;

// Load an identity matrix
void IdentityMatrixAsm( float *pfM )
{
    __asm
    {
        push        ebx
        mov			eax, pfM
        mov         ebx, 00000000h              // 0.0 in float

        mov         dword ptr[ eax + 4 ], ebx
        mov         dword ptr[ eax + 8 ], ebx
        mov         dword ptr[ eax + 12 ], ebx

        mov         dword ptr[ eax + 16 ], ebx
        mov         dword ptr[ eax + 24 ], ebx
        mov         dword ptr[ eax + 28 ], ebx

        mov         dword ptr[ eax + 32 ], ebx
        mov         dword ptr[ eax + 36 ], ebx
        mov         dword ptr[ eax + 44 ], ebx

        mov         dword ptr[ eax + 48 ], ebx
        mov         dword ptr[ eax + 52 ], ebx
        mov         dword ptr[ eax + 56 ], ebx

        mov         ebx, 3f800000h              // 1.0 in float
        mov         dword ptr[ eax ], ebx
        mov         dword ptr[ eax + 20 ], ebx
        mov         dword ptr[ eax + 40 ], ebx
        mov         dword ptr[ eax + 60 ], ebx

        pop ebx
    }
}


// loop unrolled matrix multiplication
// Basically, I load a single row from Matrix 1, push it onto the floating point stack
// and then get the vector product each column from Matrix 2 to get a row of results.
// The assembly is no longer "straightforward", I optimized by grouping reads so
// the their use is delayed as long as possible. (prevents lost cycles)
// Since the initial row values are cached in the fp stack, it IS legal to do:
// MulMatrixAsm( M1, M1, M2 ); for matrix stack operations
// MulMatrixAsm( M2, M1, M2 ); is not possible though.
void MultMatrixAsm( float *pfResults, const float *pfM1, const float *pfM2 )
{
    // assertions dont show up in release builds, so it incurs no performance loss.
    assert( pfResults != pfM2 );        

    __asm
    {
        push        ebx
        push        edx

        mov         eax, pfM1
        mov         ebx, pfM2
        mov         edx, pfResults

        
        fld         dword ptr[ eax ]
        fld         dword ptr[ eax + 4 ]
        fld         dword ptr[ eax + 8 ]
        fld         dword ptr[ eax + 12 ]
        
        fld         dword ptr[ ebx ]
        fmul        st(0), st(4)
        fld         dword ptr[ ebx + 4 ]
        fmul        st(0), st(5)
        fld         dword ptr[ ebx + 16 ]         
        fld         dword ptr[ ebx + 20 ]
        fmul        st(0), st(6)
        faddp       st(2), st(0)
        fmul        st(0), st(5)
        faddp       st(2), st(0)
        fld         dword ptr[ ebx + 32 ]         
        fld         dword ptr[ ebx + 36 ]
        fmul        st(0), st(5)
        faddp       st(2), st(0)
        fmul        st(0), st(4)
        faddp       st(2), st(0)
        fld         dword ptr[ ebx + 48 ]         
        fld         dword ptr[ ebx + 52 ]
        fmul        st(0), st(4)
        faddp       st(2), st(0)
        fmul        st(0), st(3)
        faddp       st(2), st(0)
        fstp        dword ptr[ edx + 4 ]
        fstp        dword ptr[ edx ]


        fld         dword ptr[ ebx + 12 ]
        fmul        st(0), st(4)
        fld         dword ptr[ ebx + 8 ]
        fmulp       st(5), st(0)
        fld         dword ptr[ ebx + 24 ]
        fld         dword ptr[ ebx + 28 ]
        fmul        st(0), st(5)
        faddp       st(2), st(0)
        fmulp       st(4), st(0)
        fld         dword ptr[ ebx + 40 ]
        fld         dword ptr[ ebx + 44 ]
        fmul        st(0), st(4)
        faddp       st(2), st(0)
        fmulp       st(3), st(0)
        fld         dword ptr[ ebx + 56 ]
        fld         dword ptr[ ebx + 60 ]         
        fmul        st(0), st(3)
        faddp       st(2), st(0)
        fmulp       st(2), st(0)
        fstp        dword ptr[ edx + 12 ]
        faddp       st(1), st(0)
        faddp       st(1), st(0)
        faddp       st(1), st(0)
        fstp        dword ptr[ edx + 8 ]





        fld         dword ptr[ eax + 16 ]
        fld         dword ptr[ eax + 20 ]
        fld         dword ptr[ eax + 24 ]
        fld         dword ptr[ eax + 28 ]
        
        fld         dword ptr[ ebx ]
        fmul        st(0), st(4)
        fld         dword ptr[ ebx + 4 ]
        fmul        st(0), st(5)
        fld         dword ptr[ ebx + 16 ]         
        fld         dword ptr[ ebx + 20 ]
        fmul        st(0), st(6)
        faddp       st(2), st(0)
        fmul        st(0), st(5)
        faddp       st(2), st(0)
        fld         dword ptr[ ebx + 32 ]         
        fld         dword ptr[ ebx + 36 ]
        fmul        st(0), st(5)
        faddp       st(2), st(0)
        fmul        st(0), st(4)
        faddp       st(2), st(0)
        fld         dword ptr[ ebx + 48 ]         
        fld         dword ptr[ ebx + 52 ]
        fmul        st(0), st(4)
        faddp       st(2), st(0)
        fmul        st(0), st(3)
        faddp       st(2), st(0)
        fstp        dword ptr[ edx + 20 ]
        fstp        dword ptr[ edx + 16 ]


        fld         dword ptr[ ebx + 12 ]
        fmul        st(0), st(4)
        fld         dword ptr[ ebx + 8 ]
        fmulp       st(5), st(0)
        fld         dword ptr[ ebx + 24 ]
        fld         dword ptr[ ebx + 28 ]
        fmul        st(0), st(5)
        faddp       st(2), st(0)
        fmulp       st(4), st(0)
        fld         dword ptr[ ebx + 40 ]
        fld         dword ptr[ ebx + 44 ]
        fmul        st(0), st(4)
        faddp       st(2), st(0)
        fmulp       st(3), st(0)
        fld         dword ptr[ ebx + 56 ]
        fld         dword ptr[ ebx + 60 ]         
        fmul        st(0), st(3)
        faddp       st(2), st(0)
        fmulp       st(2), st(0)
        fstp        dword ptr[ edx + 28 ]
        faddp       st(1), st(0)
        faddp       st(1), st(0)
        faddp       st(1), st(0)
        fstp        dword ptr[ edx + 24 ]





        fld         dword ptr[ eax + 32 ]
        fld         dword ptr[ eax + 36 ]
        fld         dword ptr[ eax + 40 ]
        fld         dword ptr[ eax + 44 ]
        
        fld         dword ptr[ ebx ]
        fmul        st(0), st(4)
        fld         dword ptr[ ebx + 4 ]
        fmul        st(0), st(5)
        fld         dword ptr[ ebx + 16 ]         
        fld         dword ptr[ ebx + 20 ]
        fmul        st(0), st(6)
        faddp       st(2), st(0)
        fmul        st(0), st(5)
        faddp       st(2), st(0)
        fld         dword ptr[ ebx + 32 ]         
        fld         dword ptr[ ebx + 36 ]
        fmul        st(0), st(5)
        faddp       st(2), st(0)
        fmul        st(0), st(4)
        faddp       st(2), st(0)
        fld         dword ptr[ ebx + 48 ]         
        fld         dword ptr[ ebx + 52 ]
        fmul        st(0), st(4)
        faddp       st(2), st(0)
        fmul        st(0), st(3)
        faddp       st(2), st(0)
        fstp        dword ptr[ edx + 36 ]
        fstp        dword ptr[ edx + 32 ]


        fld         dword ptr[ ebx + 12 ]
        fmul        st(0), st(4)
        fld         dword ptr[ ebx + 8 ]
        fmulp       st(5), st(0)
        fld         dword ptr[ ebx + 24 ]
        fld         dword ptr[ ebx + 28 ]
        fmul        st(0), st(5)
        faddp       st(2), st(0)
        fmulp       st(4), st(0)
        fld         dword ptr[ ebx + 40 ]
        fld         dword ptr[ ebx + 44 ]
        fmul        st(0), st(4)
        faddp       st(2), st(0)
        fmulp       st(3), st(0)
        fld         dword ptr[ ebx + 56 ]
        fld         dword ptr[ ebx + 60 ]         
        fmul        st(0), st(3)
        faddp       st(2), st(0)
        fmulp       st(2), st(0)
        fstp        dword ptr[ edx + 44 ]
        faddp       st(1), st(0)
        faddp       st(1), st(0)
        faddp       st(1), st(0)
        fstp        dword ptr[ edx + 40 ]





        fld         dword ptr[ eax + 48 ]
        fld         dword ptr[ eax + 52 ]
        fld         dword ptr[ eax + 56 ]
        fld         dword ptr[ eax + 60 ]
        
        fld         dword ptr[ ebx ]
        fmul        st(0), st(4)
        fld         dword ptr[ ebx + 4 ]
        fmul        st(0), st(5)
        fld         dword ptr[ ebx + 16 ]         
        fld         dword ptr[ ebx + 20 ]
        fmul        st(0), st(6)
        faddp       st(2), st(0)
        fmul        st(0), st(5)
        faddp       st(2), st(0)
        fld         dword ptr[ ebx + 32 ]         
        fld         dword ptr[ ebx + 36 ]
        fmul        st(0), st(5)
        faddp       st(2), st(0)
        fmul        st(0), st(4)
        faddp       st(2), st(0)
        fld         dword ptr[ ebx + 48 ]         
        fld         dword ptr[ ebx + 52 ]
        fmul        st(0), st(4)
        faddp       st(2), st(0)
        fmul        st(0), st(3)
        faddp       st(2), st(0)
        fstp        dword ptr[ edx + 52 ]
        fstp        dword ptr[ edx + 48 ]


        fld         dword ptr[ ebx + 12 ]
        fmul        st(0), st(4)
        fld         dword ptr[ ebx + 8 ]
        fmulp       st(5), st(0)
        fld         dword ptr[ ebx + 24 ]
        fld         dword ptr[ ebx + 28 ]
        fmul        st(0), st(5)
        faddp       st(2), st(0)
        fmulp       st(4), st(0)
        fld         dword ptr[ ebx + 40 ]
        fld         dword ptr[ ebx + 44 ]
        fmul        st(0), st(4)
        faddp       st(2), st(0)
        fmulp       st(3), st(0)
        fld         dword ptr[ ebx + 56 ]
        fld         dword ptr[ ebx + 60 ]         
        fmul        st(0), st(3)
        faddp       st(2), st(0)
        fmulp       st(2), st(0)
        fstp        dword ptr[ edx + 60 ]
        faddp       st(1), st(0)
        faddp       st(1), st(0)
        faddp       st(1), st(0)
        fstp        dword ptr[ edx + 56 ]



        pop         edx
        pop         ebx
    }
}


/*
// This is the loop version of the above...
__forceinline void MultMatrixAsm( float *pfResults, const float *pfM1, const float *pfM2 )
{
    // assertions dont show up in release builds, so it incurs no performance loss.
    assert( Results != M2 );        

    __asm
    {
        push        ebx
        push        ecx
        push        edx

        mov         eax, pfM1
        mov         ebx, pfM2
        mov         edx, pfResults
        mov         ecx, 48



row:
        fld         dword ptr[ eax + ecx ]
        fld         dword ptr[ eax + ecx + 4 ]
        fld         dword ptr[ eax + ecx + 8 ]
        fld         dword ptr[ eax + ecx + 12 ]
        

        fld         dword ptr[ ebx ]
        fmul        st(0), st(4)
        fld         dword ptr[ ebx + 4 ]
        fmul        st(0), st(5)
        fld         dword ptr[ ebx + 16 ]         
        fld         dword ptr[ ebx + 20 ]
        fmul        st(0), st(6)
        faddp       st(2), st(0)
        fmul        st(0), st(5)
        faddp       st(2), st(0)
        fld         dword ptr[ ebx + 32 ]         
        fld         dword ptr[ ebx + 36 ]
        fmul        st(0), st(5)
        faddp       st(2), st(0)
        fmul        st(0), st(4)
        faddp       st(2), st(0)
        fld         dword ptr[ ebx + 48 ]         
        fld         dword ptr[ ebx + 52 ]
        fmul        st(0), st(4)
        faddp       st(2), st(0)
        fmul        st(0), st(3)
        faddp       st(2), st(0)
        fstp        dword ptr[ edx + ecx + 4 ]
        fstp        dword ptr[ edx + ecx ]

        
        fld         dword ptr[ ebx + 12 ]
        fmul        st(0), st(4)
        fld         dword ptr[ ebx + 8 ]
        fmulp       st(5), st(0)
        fld         dword ptr[ ebx + 24 ]
        fld         dword ptr[ ebx + 28 ]
        fmul        st(0), st(5)
        faddp       st(2), st(0)
        fmulp       st(4), st(0)
        fld         dword ptr[ ebx + 40 ]
        fld         dword ptr[ ebx + 44 ]
        fmul        st(0), st(4)
        faddp       st(2), st(0)
        fmulp       st(3), st(0)
        fld         dword ptr[ ebx + 56 ]
        fld         dword ptr[ ebx + 60 ]         
        fmul        st(0), st(3)
        faddp       st(2), st(0)
        fmulp       st(2), st(0)
        fstp        dword ptr[ edx + ecx + 12 ]
        faddp       st(1), st(0)
        faddp       st(1), st(0)
        faddp       st(1), st(0)
        fstp        dword ptr[ edx + ecx + 8 ]

        sub         ecx, 16
        jge         row
                    
        pop         edx
        pop         ecx
        pop         ebx
    }
}
*/


// C Version.. VC7 can loop unroll this, but its not that much faster still..
// M1, M2, and Results must be different pointers, a version could be written
void CMultMatrix( float *pfResults, float *pfM1, float *pfM2 )
{
    // assertions dont show up in release builds, so it incurs no performance loss.
    assert( pfResults != pfM1 );
    assert( pfResults != pfM2 );

    for( int i = 0; i < 16; i+=4 )
    {
        for ( int j = 0; j < 4; j++ )
            pfResults[ i + j ] = pfM1[ i ] * pfM2[ j ] + 
                                 pfM1[ i + 1 ] * pfM2[ j + 4 ] + 
                                 pfM1[ i + 2 ] * pfM2[ j + 8 ] + 
                                 pfM1[ i + 3 ] * pfM2[ j + 12 ];
    }
}

void PrintMatrix( float * M )
{
    cout << endl;
    for ( int i = 0; i < 16; i+=4 )
    {
        for ( int j = 0; j < 4; j++ )
            cout << M[ i + j ] << "\t";
        cout << endl;
    }
}

int __cdecl main(int argc, char* argv[])
{
	float M1[16];
	float M2[16];
    float Results[16];
    int i;

    for( i = 0; i < 16; i++ )
    {
        M1[ i ] = 16.0f - i;
        M2[ i ] = i;
    }
    LARGE_INTEGER CTestStart, CTestFinish;
    LARGE_INTEGER D3DTestStart, D3DTestFinish;
    LARGE_INTEGER AsmTestStart, AsmTestFinish;

    // Run CMultMatrix test
    memset( Results, 0, sizeof( float ) * 16 );
    QueryPerformanceCounter( &CTestStart );
    for ( i = 0; i < 5000000; i++ )
        CMultMatrix( Results, M1, M2 );
    QueryPerformanceCounter( &CTestFinish );
    PrintMatrix( Results );
    LARGE_INTEGER CTestTime;
    CTestTime.QuadPart = CTestFinish.QuadPart - CTestStart.QuadPart;
    // number too big to print
    assert( CTestTime.HighPart == 0 );
    cout << "C  :\t" << CTestTime.LowPart << endl;

    // Run D3DXMatrixMultiply test
    D3DXMATRIX D3DResults;
    D3DXMATRIX D3DM1(M1);
    D3DXMATRIX D3DM2(M2);
    memset( D3DResults, 0, sizeof( float ) * 16 );
    QueryPerformanceCounter( &D3DTestStart );
    for ( i = 0; i < 5000000; i++ )
        D3DXMatrixMultiply( &D3DResults, &D3DM1, &D3DM2 );
    QueryPerformanceCounter( &D3DTestFinish );
    PrintMatrix( D3DResults );
    LARGE_INTEGER D3DTestTime;
    D3DTestTime.QuadPart = D3DTestFinish.QuadPart - D3DTestStart.QuadPart;
    // number too big to print
    assert( D3DTestTime.HighPart == 0 );
    cout << "D3D:\t" << D3DTestTime.LowPart << endl;

    // Run MultMatrixAsm test
    memset( Results, 0, sizeof( float ) * 16 );
    QueryPerformanceCounter( &AsmTestStart );
    for ( i = 0; i < 5000000; i++ )
        MultMatrixAsm( Results, M1, M2 );
    QueryPerformanceCounter( &AsmTestFinish );
    PrintMatrix( Results );
    LARGE_INTEGER AsmTestTime;
    AsmTestTime.QuadPart = AsmTestFinish.QuadPart - AsmTestStart.QuadPart;
    // number too big to print
    assert( AsmTestTime.HighPart == 0 );
    cout << "Asm:\t" << AsmTestTime.LowPart << endl;

    getchar();

	return 0;
}

Currently browsing [multmatrixasm.zip] (4,499 bytes) - [StdAfx.cpp] - (294 bytes)

// stdafx.cpp : source file that includes just the standard includes
//	matmult.pch will be the pre-compiled header
//	stdafx.obj will contain the pre-compiled type information

#include "stdafx.h"

// TODO: reference any additional headers you need in STDAFX.H
// and not in this file

Currently browsing [multmatrixasm.zip] (4,499 bytes) - [StdAfx.h] - (769 bytes)

// stdafx.h : include file for standard system include files,
//  or project specific include files that are used frequently, but
//      are changed infrequently
//

#if !defined(AFX_STDAFX_H__8F4EF30C_35AE_46C4_B626_F3D0050F8467__INCLUDED_)
#define AFX_STDAFX_H__8F4EF30C_35AE_46C4_B626_F3D0050F8467__INCLUDED_

#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000

#define WIN32_LEAN_AND_MEAN		// Exclude rarely-used stuff from Windows headers

#include <stdio.h>

// TODO: reference additional headers your program requires here

//{{AFX_INSERT_LOCATION}}
// Microsoft Visual C++ will insert additional declarations immediately before the previous line.

#endif // !defined(AFX_STDAFX_H__8F4EF30C_35AE_46C4_B626_F3D0050F8467__INCLUDED_)

The zip file viewer built into the Developer Toolbox made use of the zlib library, as well as the zlibdll source additions.

 

Copyright 1999-2008 (C) FLIPCODE.COM and/or the original content author(s). All rights reserved.
Please read our Terms, Conditions, and Privacy information.