Re: [Bug-gnubg] Vectorizing 3rd step

[Øystein Johansen]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

Øystein Johansen wrote:
| Here's a patch for i386 / GCC vectorizing of the inner loops of
| Evaluate(). I see some improvement, but I believe this can be
| improved even further.
|
| Some comments: I believe having a integer counter in the loop slows
| it down. Can I exit the loop in an other way?
|
| I initialize a vector for scaling in the second loop. I believe this
| can be made simpler. Any suggestions?
|
| Please comment on this two issues.

Today I've aligned the arrays in memory. That gained some speed as well.
I still believe it can be improved a bit further. My GCC 3.4.2 builds
gives me now 25200 eval/sec on this computer. I used to have about 20500
eval/sec. I'm still missing a bit 30000 I get with Ingos ICC build, but
I believe the way to go is with intrinsics.

Aligned patch attached. There is also a mm_malloc.h file you need if
you're not using a gcc snapshot.

- -Øystein

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.2.4 (MingW32)
Comment: Using GnuPG with Thunderbird - http://enigmail.mozdev.org

iD8DBQFCZC4G6kDTFPhwyqYRAhI8AKCLmdmXX3gJivdbutrmcJejaGaVZwCeKEif
rfRCRCxPuE0EAyuuZLQ1NhU=
=fNcM
-----END PGP SIGNATURE-----

Index: neuralnet.c
===================================================================
RCS file: /cvsroot/gnubg/gnubg/lib/neuralnet.c,v
retrieving revision 1.23
diff -u -r1.23 neuralnet.c
--- neuralnet.c 25 Feb 2005 11:34:24 -0000      1.23
+++ neuralnet.c 18 Apr 2005 22:04:28 -0000
@@ -34,6 +34,8 @@
 #include <time.h>
 #include <assert.h>
 
+#include <mm_malloc.h>
+
 #define SIGMOID_BAUR 0
 #define SIGMOID_JTH 0
 
@@ -342,33 +344,33 @@
     pnn->nTrained = 0;
     pnn->fDirect = FALSE;
    
-    if( !( pnn->arHiddenWeight = malloc( cHidden * cInput *
-                                        sizeof( float ) ) ) )
+    if( !( pnn->arHiddenWeight = _mm_malloc( cHidden * cInput *
+                                        sizeof( float ), 16 ) ) )
        return -1;
 
-    if( !( pnn->arOutputWeight = malloc( cOutput * cHidden *
-                                        sizeof( float ) ) ) ) {
-       free( pnn->arHiddenWeight );
+    if( !( pnn->arOutputWeight = _mm_malloc( cOutput * cHidden *
+                                        sizeof( float ), 16 ) ) ) {
+       _mm_free( pnn->arHiddenWeight );
        return -1;
     }
     
-    if( !( pnn->arHiddenThreshold = malloc( cHidden * sizeof( float ) ) ) ) {
-       free( pnn->arOutputWeight );
-       free( pnn->arHiddenWeight );
+    if( !( pnn->arHiddenThreshold = _mm_malloc( cHidden * sizeof( float ), 16 
) ) ) {
+       _mm_free( pnn->arOutputWeight );
+       _mm_free( pnn->arHiddenWeight );
        return -1;
     }
           
-    if( !( pnn->arOutputThreshold = malloc( cOutput * sizeof( float ) ) ) ) {
-       free( pnn->arHiddenThreshold );
-       free( pnn->arOutputWeight );
-       free( pnn->arHiddenWeight );
+    if( !( pnn->arOutputThreshold = _mm_malloc( cOutput * sizeof( float ), 16 
) ) ) {
+       _mm_free( pnn->arHiddenThreshold );
+       _mm_free( pnn->arOutputWeight );
+       _mm_free( pnn->arHiddenWeight );
        return -1;
     }
 
     CheckRC();
 
-    pnn->savedBase = malloc( cHidden * sizeof( float ) ); 
-    pnn->savedIBase = malloc( cInput * sizeof( float ) ); 
+    pnn->savedBase = _mm_malloc( cHidden * sizeof( float ), 16 ); 
+    pnn->savedIBase = _mm_malloc( cInput * sizeof( float ), 16 ); 
  
     for( i = cHidden * cInput, pf = pnn->arHiddenWeight; i; i-- )
        *pf++ = ( (int) ( irand( &rc ) & 0xFFFF ) - 0x8000 ) / 131072.0;
@@ -421,8 +423,8 @@
    pnn->arOutputThreshold = fp;
    fp += pnn->cOutput;
 
-   pnn->savedBase = malloc( pnn->cHidden * sizeof( float ) ); 
-   pnn->savedIBase = malloc( pnn->cInput * sizeof( float ) ); 
+   pnn->savedBase = _mm_malloc( pnn->cHidden * sizeof( float ), 16 ); 
+   pnn->savedIBase = _mm_malloc( pnn->cInput * sizeof( float ), 16 ); 
 
    return fp;
 }
@@ -432,18 +434,33 @@
 NeuralNetDestroy( neuralnet *pnn )
 {
   if( !pnn->fDirect ) {
-    free( pnn->arHiddenWeight ); pnn->arHiddenWeight = 0;
-    free( pnn->arOutputWeight ); pnn->arOutputWeight = 0;
-    free( pnn->arHiddenThreshold ); pnn->arHiddenThreshold = 0;
-    free( pnn->arOutputThreshold ); pnn->arOutputThreshold = 0;
+    _mm_free( pnn->arHiddenWeight ); pnn->arHiddenWeight = 0;
+    _mm_free( pnn->arOutputWeight ); pnn->arOutputWeight = 0;
+    _mm_free( pnn->arHiddenThreshold ); pnn->arHiddenThreshold = 0;
+    _mm_free( pnn->arOutputThreshold ); pnn->arOutputThreshold = 0;
   }
 
-  free(pnn->savedBase); pnn->savedBase = 0;
-  free(pnn->savedIBase); pnn->savedIBase = 0;
+  _mm_free(pnn->savedBase); pnn->savedBase = 0;
+  _mm_free(pnn->savedIBase); pnn->savedIBase = 0;
   
   return 0;
 }
 
+typedef float v4sf __attribute__ ((vector_size(16)));
+
+typedef union _vec4f {
+  v4sf v;
+  float f[4];
+} vec4f;
+
+#if DEBUG
+void
+printvec( v4sf vec ){
+  float *pFloat = ( float *) &vec;
+  printf("%f, %f, %f, %f\n", pFloat[0], pFloat[1], pFloat[2], pFloat[3]);
+}
+#endif
+
 static int Evaluate( neuralnet *pnn, float arInput[], float ar[],
                         float arOutput[], float *saveAr ) {
 
@@ -452,6 +469,9 @@
 #else
     int i, j;
     float *prWeight;
+    v4sf sum, vec0, vec1, vec3;
+    
+    assert(pnn->cHidden == 128);
 
     /* Calculate activity at hidden nodes */
     for( i = 0; i < pnn->cHidden; i++ )
@@ -466,11 +486,28 @@
            float *pr = ar;
 
            if( ari == 1.0f )
-               for( j = pnn->cHidden; j; j-- )
-                   *pr++ += *prWeight++;
-           else
-               for( j = pnn->cHidden; j; j-- )
-                   *pr++ += *prWeight++ * ari;
+               for( j = 32; j; j--, pr += 4, prWeight += 4 ){
+                   vec0 = __builtin_ia32_loadaps(pr);  
+                  vec1 = __builtin_ia32_loadaps(prWeight); 
+                   sum = __builtin_ia32_addps(vec0, vec1);
+                   __builtin_ia32_storeaps (pr, sum);
+               }
+//                 *pr++ += *prWeight++;
+           else {
+               float scale[4];
+               v4sf scalevector;
+                scale[0] = scale[1] = scale[2] = scale[3] = ari;
+               scalevector = __builtin_ia32_loadaps(scale);
+               for( j = 32; j; j--, pr += 4, prWeight += 4 ){
+                   vec0 = __builtin_ia32_loadaps(pr);  
+                  vec1 = __builtin_ia32_loadaps(prWeight); 
+                  vec3 = __builtin_ia32_mulps(vec1, scalevector);
+                   sum = __builtin_ia32_addps(vec0, vec3);
+                   __builtin_ia32_storeaps (pr, sum);
+               }
+//             for( j = pnn->cHidden; j; j-- )
+//                 *pr++ += *prWeight++ * ari;
+           }
        } else
            prWeight += pnn->cHidden;
     }
@@ -484,14 +521,22 @@
 
     /* Calculate activity at output nodes */
     prWeight = pnn->arOutputWeight;
-
+    
     for( i = 0; i < pnn->cOutput; i++ ) {
-       float r = pnn->arOutputThreshold[ i ];
-       
-       for( j = 0; j < pnn->cHidden; j++ )
-           r += ar[ j ] * *prWeight++;
-
-       arOutput[ i ] = sigmoid( -pnn->rBetaOutput * r );
+       float r = pnn->arOutputThreshold[ i ];
+       float *pr = ar;
+       vec4f sum;
+       v4sf vec0, vec1, vec3;
+       sum.v = __builtin_ia32_xorps(sum.v, sum.v);
+       for( j = 32; j ; j--, prWeight += 4, pr += 4 ){
+         vec0 = __builtin_ia32_loadaps(pr);       /* Four floats into vec0 */
+         vec1 = __builtin_ia32_loadaps(prWeight); /* Four weights into vect1 
*/ 
+         vec3 = __builtin_ia32_mulps(vec0, vec1); /* Multiply */
+         sum.v = __builtin_ia32_addps(sum.v, vec3); /* Add */
+       }
+       
+       r += sum.f[0] + sum.f[1] + sum.f[2] + sum.f[3]; 
+       arOutput[ i ] = sigmoid( -pnn->rBetaOutput * r );
     }
 
     return 0;
@@ -879,7 +924,7 @@
     }
     
     if( cHidden != pnn->cHidden || cInput != pnn->cInput ) {
-       if( !( pr = malloc( cHidden * cInput * sizeof( float ) ) ) )
+       if( !( pr = _mm_malloc( cHidden * cInput * sizeof( float ), 16 ) ) )
            return -1;
 
        prNew = pr;
@@ -895,7 +940,7 @@
                else
                    *prNew++ = pnn->arHiddenWeight[ i * pnn->cHidden + j ];
                    
-       free( pnn->arHiddenWeight );
+       _mm_free( pnn->arHiddenWeight );
 
        pnn->arHiddenWeight = pr;
     }
@@ -911,7 +956,7 @@
     }
     
     if( cOutput != pnn->cOutput || cHidden != pnn->cHidden ) {
-       if( !( pr = malloc( cOutput * cHidden * sizeof( float ) ) ) )
+       if( !( pr = _mm_malloc( cOutput * cHidden * sizeof( float ), 16 ) ) )
            return -1;
 
        prNew = pr;
@@ -927,7 +972,7 @@
                else
                    *prNew++ = pnn->arOutputWeight[ i * pnn->cOutput + j ];
 
-       free( pnn->arOutputWeight );
+       _mm_free( pnn->arOutputWeight );
 
        pnn->arOutputWeight = pr;
     }

/* Copyright (C) 2004 Free Software Foundation, Inc.

   This file is part of GCC.

   GCC is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2, or (at your option)
   any later version.

   GCC is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with GCC; see the file COPYING.  If not, write to
   the Free Software Foundation, 59 Temple Place - Suite 330,
   Boston, MA 02111-1307, USA.  */

/* As a special exception, if you include this header file into source
   files compiled by GCC, this header file does not by itself cause
   the resulting executable to be covered by the GNU General Public
   License.  This exception does not however invalidate any other
   reasons why the executable file might be covered by the GNU General
   Public License.  */

#ifndef _MM_MALLOC_H_INCLUDED
#define _MM_MALLOC_H_INCLUDED

#include <stdlib.h>
#include <errno.h>

static __inline__ void*
_mm_malloc (size_t size, size_t align)
{
  void * malloc_ptr;
  void * aligned_ptr;

  /* Error if align is not a power of two.  */
  if (align & (align - 1))
    {
      errno = EINVAL;
      return ((void*) 0);
    }

  if (size == 0)
    return ((void *) 0);

 /* Assume malloc'd pointer is aligned at least to sizeof (void*).
    If necessary, add another sizeof (void*) to store the value
    returned by malloc. Effectively this enforces a minimum alignment
    of sizeof double. */
    if (align < 2 * sizeof (void *))
      align = 2 * sizeof (void *);

  malloc_ptr = malloc (size + align);
  if (!malloc_ptr)
    return ((void *) 0);

  /* Align  We have at least sizeof (void *) space below malloc'd ptr. */
  aligned_ptr = (void *) (((size_t) malloc_ptr + align)
                          & ~((size_t) (align) - 1));

  /* Store the original pointer just before p.  */
  ((void **) aligned_ptr) [-1] = malloc_ptr;

  return aligned_ptr;
}

static __inline__ void
_mm_free (void * aligned_ptr)
{
  if (aligned_ptr)
    free (((void **) aligned_ptr) [-1]);
}

#endif /* _MM_MALLOC_H_INCLUDED */