mini-SSE-L1-BLAS library:
A fast library for SDOT,DDOT,SAXPY,DAXPY operations on x86 processor

The mini-SSE-L1-BLAS Library is providing the basic Level 1 BLAS operations on vectors. In particular, the four functions: SDOT, DDOT, SAXPY, DAXPY were deeply hand-optimized. The library is composed of only 2 small files (a ".cpp" and a ".h" file). Inclusion of the library inside your own projects is easy. It is cross-platform. It compiles transparently on both Visual Studio .NET and GCC. A FORTRAN interface is also available. Optimization includes inline SSE and SSE2 assembly code inserted inside the C/C++ code, loop unrolling. Assembly code instructions are ordered to increase parallel execution of instructions, to ease branch prediction, to reduce dependency links between two close instructions. Software Prefetch instructions were very roughly investigated but were not included in the library. The functions are also implemented as C++ macros to remove the function calls overhead.

A complete User's guide for the FORTRAN interface and for the C++ interface is given respectively in section 1.3 and section 1.4. If you are in a hurry, you can skip all the other sections and only read the user's guides (sections 1.3 and 1.4)..

• As a general rule of thumb, when optimizing code you should never try to guess the best strategy based on some empirical reasoning. Always validate your choices using experimental benchmarks results only.
  
  
• SSE and SSE2 are a set of assembly-level instructions that allows to compute 4 products of float's or 2 products of double's in, essentially, the same time as 1 product using the classical set of assembly-level instructions. This normally leads to a reduction of the computation time for a SDOT or SAXPY operation to 25% of a non-optimized classical implementation. For a DDOT or DAXPY operation the computing time should be reduced to 50% . Unfortunately, the real bottleneck for high performances is memory access. This means that we cannot usually obtain such good performance but still we can get close enough.
  
  
• Let's consider this small C code:

  int i=n;
  double s=0.0;
  while (i--) s+=a[i]*b[i];
  

  This code computes the dot-product of 2 vectors of length n: a and b. To reduce computation time you can unroll (8 times) the main loop:

  int i=n>>3, j=0, k=n&7;
  double s=0.0;
  while (i--)
  {
      s+=a[j  ]*b[j  ];
      s+=a[j+1]*b[j+1];
      s+=a[j+2]*b[j+2];
      s+=a[j+3]*b[j+3];
      s+=a[j+4]*b[j+4];
      s+=a[j+5]*b[j+5];
      s+=a[j+6]*b[j+6];
      s+=a[j+7]*b[j+7];
      j+=8;
  }
  switch(k)
  {
      case 7: s+=a[j]*b[j]; j++;
      case 6: s+=a[j]*b[j]; j++;
      case 5: s+=a[j]*b[j]; j++;
      case 4: s+=a[j]*b[j]; j++;
      case 3: s+=a[j]*b[j]; j++;
      case 2: s+=a[j]*b[j]; j++;
      case 1: s+=a[j]*b[j]; j++;
  }
  

  The second code is faster because it performs the end-of-loop test "if (i==0)" once every height multiplies.
  
  
• Let's consider this small C code (we assume that and ):

  int i=0;
  register double at,bt,st,s=0.0;
  while (i<n)
  {
      at=a[i  ]; bt=b[i  ]; st=at*bt; s+=st;
      at=a[i+1]; bt=b[i+1]; st=at*bt; s+=st;
      i+=2;
  }
  

  To reduce computation time you can re-order the instruction inside the main loop:

  int i=0;
  register double at,bt,st,s=0.0;
  n-=2;
  bt=b[i];
  while (i<n)
  {
      at=a[i  ]; st=at*bt; bt=b[i+1]; s+=st;
      at=a[i+1]; st=at*bt; bt=b[i+3]; s+=st;
      i+=2;
  }
  s+= a[i  ] * bt    +
      a[i+1] * b[i+1];
  

  The second code is faster because it allows parallel execution of instructions. While the arithmetic-unit of the processor is busy computing "st=at*bt" the memory-access-unit of the processor can already performs "bt=b[i+1]". This is called "parallel execution" because the two units are working in parallel. Unfortunately, in our example, the next instruction (after "bt=b[i+1]") is "s+=st". To execute this instruction we must wait until the multiplication "st=at*bt" is finished, loosing some precious time. In other words, there is a "dependency link" between the instructions "st=at*bt" and "s+=st". We should re-order the instructions to remove those "dependency link" to increase parallel execution.
  
  
• Let's consider again the unrolled dot product:

  int i=n>>3, j=0, k=n&7;
  double s=0.0;
  while (i--)
  {
      s+=a[j  ]*b[j  ];
      s+=a[j+1]*b[j+1];
      s+=a[j+2]*b[j+2];
      s+=a[j+3]*b[j+3];
      s+=a[j+4]*b[j+4];
      s+=a[j+5]*b[j+5];
      s+=a[j+6]*b[j+6];
      s+=a[j+7]*b[j+7];
      j+=8;
  }
  switch(k)
  {
      case 7: s+=a[j]*b[j]; j++;
      case 6: s+=a[j]*b[j]; j++;
      case 5: s+=a[j]*b[j]; j++;
      case 4: s+=a[j]*b[j]; j++;
      case 3: s+=a[j]*b[j]; j++;
      case 2: s+=a[j]*b[j]; j++;
      case 1: s+=a[j]*b[j]; j++;
  }
  

  There are different ways to implement this "switch instruction" in assembler:
  
  1. CHOICE_1:

                     k=n&7;
                     goto JUMP_TABLE[k];
     JUMP_TABLE[7]:  s+=a[j]*b[j]; j++;
     JUMP_TABLE[6]:  s+=a[j]*b[j]; j++;
     JUMP_TABLE[5]:  s+=a[j]*b[j]; j++;
     JUMP_TABLE[4]:  s+=a[j]*b[j]; j++;
     JUMP_TABLE[3]:  s+=a[j]*b[j]; j++;
     JUMP_TABLE[2]:  s+=a[j]*b[j]; j++;
     JUMP_TABLE[1]:  s+=a[j]*b[j]; j++;
     JUMP_TABLE[0]:
     

     Note that the C code used here will not actually compile: It's just to have a general idea of what is possible in assembler.
     
  2. CHOICE_2:

     k=n&7;
     if (k)
     {
         s+=a[j]*b[j]; j++; k--;
         if (k)
         {
             s+=a[j]*b[j]; j++; k--;
             if (k)
             {
                 s+=a[j]*b[j]; j++; k--;
                 if (k)
                 {
                     s+=a[j]*b[j]; j++; k--;
                     if (k)
                     {
                         s+=a[j]*b[j]; j++; k--;
                         if (k)
                         {
                             s+=a[j]*b[j]; j++; k--;
                             if (k)
                             {
                                 s+=a[j]*b[j]; j++;
                             }
                         }
                     }
                 }
             }
         }
     }
     

  3. CHOICE_3:

     k=n&7;
     if (k)
     {
         s+=a[j]*b[j];
         if (k>3)
         {
             s+=a[j+1]*b[j+1];
             s+=a[j+2]*b[j+2];
             s+=a[j+3]*b[j+3];
             if (k>5)
             {
                 s+=a[j+4]*b[j+4];
                 s+=a[j+5]*b[j+5];
                 if (k==7) s+=a[j+6]*b[j+6];
             } else
             {
                 if (k==5) s+=a[j+4]*b[j+4];
             }
         } else
         {
             if (k==3)
             {
                 s+=a[j+1]*b[j+1];
                 s+=a[j+2]*b[j+2];
             } else
             {
                 if (k==2)
                 {
                     s+=a[j+1]*b[j+1];
                 }
             }
         }
         j+=k;
     }
     

  4. CHOICE_4:

                     k=n&7;
                     double ss=0.0;
                     goto JUMP_TABLE[k];
     JUMP_TABLE[7]:  s +=a[j+6]*b[j+6];
     JUMP_TABLE[6]:  ss+=a[j+5]*b[j+5];
     JUMP_TABLE[5]:  s +=a[j+4]*b[j+4];
     JUMP_TABLE[4]:  ss+=a[j+3]*b[j+3];
     JUMP_TABLE[3]:  s +=a[j+2]*b[j+2];
     JUMP_TABLE[2]:  ss+=a[j+1]*b[j+1];
     JUMP_TABLE[1]:  s +=a[j  ]*b[j];
                     j+=k;
                     s+=ss;
     JUMP_TABLE[0]:
     

  All the compilers will use either CHOICE_1 or CHOICE_2. In particular, CHOICE_1 *seems* very promising because the processor has no test (no "if" instruction) to perform.
  
  As previously mentioned, the processor is always trying to execute instructions in parallel. This means that while it is processing the current instruction it's already preparing and executing instruction "ahead". When the processor sees a GOTO or JUMP it should be able to guess where it will go using a technique called "branch prediction". A correct guess allows to prepare "ahead" the execution of instructions. A wrong guess will cause a major stall in the pipeline of instructions. Recovering from this pipeline stall is very time-consuming.
  
  CHOICE_3 is superior to CHOICE_2 because it involves less tests. CHOICE_4 is superior to CHOICE_1 because there are less dependency between the instructions. CHOICE_3 and CHOICE_4 were both implemented. Benchmark results indicates that CHOICE_3 is the fastest. CHOICE_3 was thus thoroughly used inside the library. It appears that the only jump performed in CHOICE_4 is harder to predict than all the small jumps performed in CHOICE_3.
  
  
• Variables stored in memory are slow to access. This is why there are two levels of cache on x86 processors. The following table summarizes the time required to read data depending on their location:
  
  

  Location of Data	Read time
  L1	< 3 nS
  L2	< 10 nS
  RAM	<100 nS assuming no page page misses plus possible delays to write back a dirty cache line
  Disk	10+ milliseconds
  Network Disk	100 mS to tens of seconds

  
  The objective of the Prefetch instructions is to load in advance some data in the L1 or L2 cache to decrease the time needed later to read these data. The Prefetch instruction takes as parameter a memory address. The data at this address are copied into the cache. This parameter is difficult to tune. Some preliminary results show that the usage of the Prefetch instruction inside the SDOT, DDOT, SAXPY, DAXPY functions did not lead to any speed increase (on the constrary). I think that more tuning of the Prefetch parameters should be performed. Prefetch instructions are thus currently not used.

Experimental results

Annexe

double  ddot_C(int n, double *x,double *y)
    { double s=0.0; while (n--)     { s+=*(x++) * *(y++); }; return s; }

float   sdot_C(int n, float  *x,float  *y)
    { float  s=0.0; while (n--)     { s+=*(x++) * *(y++); }; return s; }

void   saxpy_C(int n, float   a,float  *x, float  *y)
    {               while (n--)  *(x++)+=  a    * *(y++); }

void   daxpy_C(int n, double  a,double *x, double *y)
    {               while (n--)  *(x++)+=  a    * *(y++); }

User's guide for the FORTRAN interface to the library

g++ -O -c miniSSEL1BLAS.cpp

f77 -o <name_of_the_executable> <list_of_fortran_.f_files> miniSSEL1BLAS.o

    integer          function   hasSSE()                 // out=1 if SSE  supported, out=0 otherwise

    integer          function   hasSSE2()                // out=1 if SSE2 supported, out=0 otherwise


                     subroutine scopy    (n,  sx,sy)     // y[i]=x[i]        for_all_i
                     subroutine dcopy    (n,  dx,dy)

                     subroutine saxpy    (n,sa,sx,sy)    // y[i]=a*x[i]+y[i] for_all_i
                     subroutine saxpyUA  (n,sa,sx,sy)
                     subroutine daxpy    (n,da,dx,dy)
                     subroutine daxpyUA  (n,da,dx,dy)

                     subroutine sscale   (n,sa,sx)       // x[i]=a*x[i]      for_all_i
                     subroutine dscale   (n,da,dx)
                     subroutine sscaleUA (n,sa,sx)
                     subroutine dscaleUA (n,da,dx)

    real             function   sdot     (n,   sx,sy)   // out = sum_over_i x[i]*y[i]
    real             function   sdotUA   (n,   sx,sy)
    double precision function   ddot     (n,   dx,dy)
    double precision function   ddotUA   (n,   dx,dy)

    real             function   ssquare  (n,   sx)      // out = sum_over_i x[i]**2
    real             function   ssquareUA(n,   sx)
    double precision function   dsquare  (n,   dx)
    double precision function   dsquareUA(n,   dx)

    real             function   snrm2    (n,   sx)      // out = sqrt(square(n,x)/n)
    real             function   snrm2UA  (n,   sx)
    double precision function   dnrm2    (n,   dx)
    double precision function   dnrm2UA  (n,   dx)

    real             function   sasum    (n,   sx)      // out = sum_over_i fabs(x[i])
    real             function   sasumUA  (n,   sx)
    double precision function   dasum    (n,   dx)
    double precision function   dasumUA  (n,   dx)

    real             function   smin     (n,   sx)      // out = min_for_all_i x[i]
    real             function   sminUA   (n,   sx)
    double precision function   dmin     (n,   dx)
    double precision function   dminUA   (n,   dx)

    real             function   smax     (n,   sx)      // out = max_for_all_i x[i]
    real             function   smaxUA   (n,   sx)
    double precision function   dmax     (n,   dx)
    double precision function   dmaxUA   (n,   dx)

                     subroutine sswap    (n,   sx,sy)   // SWAP(x[i],y[i])  for_all_i
                     subroutine dswap    (n,   dx,dy)
                     subroutine sset     (n,sa,sx)      // x[i]=a           for_all_i
                     subroutine dset     (n,da,dx)

    integer n
    real sx(*),sy(*),sa,sc,ss
    double precision dx(*),dy(*),da,dc,ds

    integer function   suiamax( n,sx)             // out = index i of the maximum of fabs(x[i])
    integer function   duiamax( n,dx)
            subroutine srot   (sx,sy,sc,ss)       // apply a 2D rotation(c,s) on all the 2D points (x[i],[i])
            subroutine drot   (dx,dy,dc,ds)
            subroutine srotg  (sa,sb,sc,ss)       // compute a 2D rotation
            subroutine drotg  (da,db,dc,ds)

        integer n
        real sx(*),sy(*)sa,sb,sc,ss
        double precision dx(*),dy(*),da,db,dc,ds

int offset,i,n
real              sx(100), sy(100), sr
double precision  dx(100), dy(100), dr


sr=sdot(n,sx,sy)
dr=ddot(n,dx,dy)

sdot(n,sx(offset),sy(offset+i*4))

ddot(n,dx(offset),dy(offset+i*2))

int offset,i,n
real s(*)
double precision d(*)

sdotUA(n,s(offset),s(offset+i*4))
ddotUA(n,d(offset),d(offset+i*2))

User's guide for the C++ interface to the library

1. a Visual Studio .NET example (.sln)
2. a GCC example (makefile based)

void *mallocAligned(unsigned long n);   // malloc aligned on 16 byte boundaries
void freeAligned(void *t);              // free the aligned allocation

int i= ...any number...
float *f=(float*)mallocAligned(...any number...),
      *fA=f+i*4;
double *d=(double*)mallocAligned(...any number...),
       *dA=d+i*2;

    float  dot     (int n,            float  *x, float  *y);  // out = sum_over_i x[i]*y[i]
    double dot     (int n,            double *x, double *y);
    void   axpy    (int n, float   a, float  *x, float  *y);  // y[i]+=a*x[i]  for_all_i
    void   axpy    (int n, double  a, double *x, double *y);
    void   vcopy   (int n,            float  *x, float  *y);  // y[i]=x[i]      for_all_i
    void   vcopy   (int n,            double *x, double *y);
    void   vscale  (int n, float   a, float  *x);             // x[i]*=a        for_all_i
    void   vscale  (int n, double  a, double *x);
    float  asum    (int n,            float  *x);             // out = sum_over_i fabs(x[i])
    double asum    (int n,            double *x);
    float  square  (int n,            float  *x);             // out = sum_over_i x[i]**2
    double square  (int n,            double *x);
    float  nrm2    (int n,            float  *x);             // out = sqrt(square(n,x)/n)
    double nrm2    (int n,            double *x);
    float  vmin    (int n,            float  *x);             // out = min_for_all_i x[i]
    double vmin    (int n,            double *x);
    float  vmax    (int n,            float  *x);             // out = max_for_all_i x[i]
    double vmax    (int n,            double *x);
    void   vswap   (int n,            float  *x, float  *y);  // SWAP(x[i],y[i])  for_all_i
    void   vswap   (int n,            double *x, double *y);
    void   vset    (int n, float   a, float  *x);             // x[i]=a           for_all_i
    void   vset    (int n, double  a, double *x);

  int  uiamax(int  n, double *x);                                 // out = index i of the maximum of fabs(x[i])
  int  uiamax(int  n, float  *x);                                    
  void rot   (int *n, float  *x, float  *y, float   c,float   s); // apply a 2D rotation(c,s) on all the 2D points (x[i],[i])
  void rot   (int *n, double *x, double *y, double  c,double  s);     
  void rotg  (        float  *a, float  *b, float  *c,float  *s); // compute a 2D rotation
  void rotg  (        double *a, double *b, double *c,double *s);

    float  dotUA   (int n,            float  *x, float  *y);
    double dotUA   (int n,            double *x, double *y);
    void   axpyUA  (int n, float   a, float  *x, float  *y);
    void   axpyUA  (int n, double  a, double *x, double *y);
    void   vscaleUA(int n, float   a, float  *x);
    void   vscaleUA(int n, double  a, double *x);
    float  asumUA  (int n,            float  *x);
    double asumUA  (int n,            double *x);
    float  squareUA(int n,            float  *x);
    double squareUA(int n,            double *x);
    float  nrm2UA  (int n,            float  *x);
    double nrm2UA  (int n,            double *x);
    float  vminUA  (int n,            float  *x);
    double vminUA  (int n,            double *x);
    float  vmaxUA  (int n,            float  *x);
    double vmaxUA  (int n,            double *x);

int i;
float  *f=(float *)malloc(...);       // f is not aligned
double *d=(double*)malloc(...);       // d is not aligned
...
sdotUA(n,f,f+i*4);       // f is not aligned
ddotUA(n,d,d+i*2);       // d is not aligned

sdot(n,f,f+i*4);         // f is not aligned
ddot(n,d,d+i*2);         // d is not aligned

unsigned int n;
    float  *xf,*yf,af,rf;  // xf,yf are correctly aligned.
    double *xd,*yd,ad,rd;  // xd,yd are correctly aligned.

       SDOT_SSE (A,n,   xf,yf,rf)    // rf = sum_over_i xf[i]*yf[i]
       DDOT_SSE2(A,n,   xd,yd,rd)    // rd = sum_over_i xd[i]*yd[i]
      SAXPY_SSE (A,n,af,xf,yf)       // yf[i]+=af*xf[i]  for_all_i
      DAXPY_SSE2(A,n,ad,xd,yd)       // yd[i]+=ad*xd[i]  for_all_i
     SSCALE_SSE (A,n,af,xf)          // xf[i]*=af        for_all_i
     DSCALE_SSE2(A,n,ad,xd)          // xd[i]*=ad        for_all_i
    SSQUARE_SSE (A,n,   xf,rf)       // rf = sum_over_i xf[i]**2
    DSQUARE_SSE2(A,n,   xd,rd)       // rd = sum_over_i xd[i]**2
      SASUM_SSE (A,n,   xd,rd)       // rf = sum_over_i fabs(xf[i])
      DASUM_SSE2(A,n,   xf,rd)       // rd = sum_over_i fabs(xd[i])
       SMIN_SSE (A,n,   xf,rf)       // rf = min_for_all_i xf[i]
       DMIN_SSE2(A,n,   xd,rd)       // rd = min_for_all_i xd[i]
       SMAX_SSE (A,n,   xf,rf)       // rf = max_for_all_i xf[i]
       DMAX_SSE2(A,n,   xd,rd)       // rd = max_for_all_i xd[i]

    float  *xfu,*yfu,af;  // xfu,yfu can be Un-Aligned.
    double *xdu,    ,ad;  // xdu     can be Un-Aligned.

    MEMSWAP_MMX (A,n,   xfu,yfu)     // SWAP(xf[i],yf[i])  for_all_i
    SMEMSET_MMX (A,n,af,xfu)         // xf[i]=af           for_all_i
    DMEMSET_MMX (A,n,ad,xdu)         // xd[i]=ad           for_all_i

...
float  *x1,*y1,r1,
       *x2,*y2,r2;
...
SDOT_SSE(A,n,x1,y1,r1);
SDOT_SSE(A,n,x2,y2,r2);
...

...
SDOT_SSE(A,n,x1,y1,r1);
SDOT_SSE(B,n,x2,y2,r2);
...

What processor can I use?

• Intel® Pentium III and above
• AMD's 3Dnow!™ Professional technology, which is supported in
  • AMD Athlon XP
  • AMD Athlon MP
  • mobile AMD Athlon 4
  • Model 7 AMD Duron™ processors.
  • AMD opteron, AMD Athlon 64 and above
  The oldest AMD processor supporting SSE is an AMD Athlon XP 1500+ at a frequency of 1333 MHz.

• Intel® Pentium 4, Intel® Pentium M
• AMD opteron, AMD Athlon 64 and above

References

• ATLAS: Automatically Tuned Linear Algebra Software:
  http://math-atlas.sourceforge.net/
  
  When you install ATLAS on your computer, it tries to guess what's the optimum loop-unrolling. It also tries different orders for the C instructions. This last option is somewhat redundant with a good compiler. The performances of ATLAS are reported to be very close to the optimum theoretical performances of the processor when using the classical set (no SSE, no SSE2) of assembly-level instructions.
  
  Recently, some x86 assembler code was included inside ATLAS. However, the coders from ATLAS don't want to bother the users with such thing like "memory alignment". In most cases, this prevents them to write SSE-optimized functions. Performance of the SDOT,DDOT,SAXPY,DAXPY functions are thus low compared to the miniSSEL1BLAS library.
  
  
• AMD Core Math Library (ACML)
  http://developer.amd.com/acml.aspx
  
  AMD is providing an optimized Matrix processing library for all its processors. This is a commercial library. Source code is not available. However free download for evaluation purposes for Windows and Linux is available. If you desire to distribute any part of this library inside a commercial product, you shall enter into a separate written agreement with AMD.
  
  
• INTEL IPP: Intel Integrated Performance Primitives
  http://www.intel.com/cd/software/products/asmo-na/eng/perflib/ipp/index.htm

  INTEL is providing an optimized Matrix processing library for its pentium processors. This is a commercial library. Source code is not available. However free download for non-commercial
  applications (and only for the Linux OS) is available.
  

• LFMat: LFMat is a templatized open source C++ matrix library with storage-compatible and assembler specializations for 3Dnow!, SSE and SSE2.
  http://lfmat.sourceforge.net/
  
  The library should give better performances than the ATLAS library since it is using SSE/SSE2 sets of instructions. It is however still in beta currently (August 2005).
  
  
• Blitz++: Scientific computing C++ library. Uses template techniques to achieve high performance.
  http://www.oonumerics.org/blitz/
  
  Blitz++ uses templates techniques to know at compile-time, for each user loop, the optimal unrolling and the optimal order of the instructions. The information needed to optimize the loops are, basically, the exact size of the vectors and the exact size of the matrices involved. These sizes are fixed and must be given at compilation time. If this is not a problem for you, Blitz++ should give even better performances than ATLAS.
  
  
• miniSSEL1BLAS library: implements only LEVEL 1 BLAS functions on x86 processor with SSE/SSE2 support.
  http://www.applied-mathematics.net/miniSSEL1BLAS/miniSSEL1BLAS.html
  
  Many assembler level optimizations have been performed to ensure that peak performances are attained. It's currently the only library that forces you to use 16-bytes-memory-aligned vector arrays. Thanks to that, it should be faster than any other implementation. And it's free! Available for direct download! Drawback: for some functions, the user MUST use 16-bytes-memory-aligned vector arrays.
  

• A good starting point: http://www.tommesani.com/Docs.html
• IA-32 Intel® Architecture Software Developer's Manual Volume 1: Basic Architecture: 24547006.pdf
• IA-32 Intel® Architecture Software Developer's Manual Volume 2: Instruction Set Reference: 24547112.pdf
• Software Optimization Guide for AMD Athlon(tm) 64 and AMD Opteron(tm) Processors: AMD_optimization.pdf
• How to optimize for the Pentium family of microprocessors by Agner Fog, Ph.D.: Optimization guide for x86 platforms (local copy here)

• James Mc Parlane's Blog: x86 Cache Prefetch:
  • http://blog.metawrap.com/blog/CategoryView,category,Assembler.aspx
• Memory Access is the Key To Speed (The table about access time for the different kind of memory is coming from this site).
  • http://www.iseran.com/Win32/CodeForSpeed/memory.html