VC++ - intrinsic optimization

VC++ does not support inline assembler when it compiles to x64 target. But it says you don't have to use inline assembler as there is intrinisics. Intrinisics has methods that match with assembler instructions.

void SSE_Pitch0( 
  const uint8_t *h_src, uint8_t *h_dst, 
  int width, int height, 
  int roiLeft, int roiTop, int roiRight, int roiBottom,
  float horPitch, float verPitch )
{
  int sizeOfXMMInWords = SSE::SizeOfDataType / 2;
  int integerPitch = (int)horPitch;

  uint16_t toFloor = (uint16_t)( (horPitch - (int)horPitch) * 0xFF);
  uint16_t toCeiling = 0xFF - toFloor;

  std::vector< uint16_t > toFloors( sizeOfXMMInWords, toFloor );
  std::vector< uint16_t > toCeilings( sizeOfXMMInWords, toCeiling );

  __m128i xmmFloor, xmmCeiling;
  xmmFloor = _mm_loadu_si128( (const __m128i*)( toFloors.front()) );
  xmmCeiling = _mm_loadu_si128( (const __m128i*)( toCeilings.front()) );

  for ( int y=roiTop; y<=roiBottom; y++ ) 
  {
    int x = roiLeft;
    const uint8_t *source = h_src + width*y+x;
    uint8_t *target = h_dst + width*y+x;

    for ( ; x<=roiRight; 
      x+=sizeOfXMMInWords, source+=sizeOfXMMInWords, target+=sizeOfXMMInWords ) 
    {
      __m128i xmmTemp, xmmLeft, xmmRight, xmmEast, xmmWest, xmmCenter;

      // _mm_cvtepu8_epi16 can't load value from memory directly unlike PMOVZXBW
      xmmTemp = _mm_loadu_si128( (const __m128i*)(source-integerPitch-1) );
      xmmEast = _mm_cvtepu8_epi16( xmmTemp );
      xmmTemp = _mm_srli_si128( xmmTemp, 1 );
      xmmWest = _mm_cvtepu8_epi16( xmmTemp );

      xmmEast = _mm_mullo_epi16( xmmEast, xmmFloor );
      xmmWest = _mm_mullo_epi16( xmmWest, xmmCeiling );
      xmmLeft = _mm_adds_epu16( xmmEast, xmmWest );
      xmmLeft = _mm_srli_epi16( xmmLeft, 1 );

      xmmTemp = _mm_loadu_si128( (const __m128i*)(source+integerPitch) );
      xmmEast = _mm_cvtepu8_epi16( xmmTemp );
      xmmTemp = _mm_srli_si128( xmmTemp, 1 );
      xmmWest = _mm_cvtepu8_epi16( xmmTemp );

      xmmEast = _mm_mullo_epi16( xmmEast, xmmCeiling );
      xmmWest = _mm_mullo_epi16( xmmWest, xmmFloor );
      xmmRight = _mm_adds_epu16( xmmEast, xmmWest );
      xmmRight = _mm_srli_epi16( xmmRight, 1 );

      xmmTemp = _mm_loadu_si128( (const __m128i*)(source) );
      xmmCenter = _mm_cvtepu8_epi16( xmmTemp );
      xmmCenter = _mm_slli_epi16( xmmCenter, 8 );
      xmmTemp = xmmCenter;

      xmmCenter = _mm_subs_epu16( xmmCenter, xmmLeft );
      xmmCenter = _mm_subs_epu16( xmmCenter, xmmRight );

      xmmLeft = _mm_adds_epu16( xmmLeft, xmmRight);
      xmmLeft = _mm_subs_epu16( xmmLeft, xmmTemp );

      xmmCenter = _mm_adds_epu16( xmmCenter, xmmLeft );
      xmmCenter = _mm_srli_epi16( xmmCenter, 8 );
      xmmTemp = _mm_xor_si128( xmmTemp, xmmTemp);
      xmmTemp = _mm_packus_epi16( xmmCenter, xmmTemp );

      _mm_storel_epi64( (__m128i*)(target), xmmTemp );
    }
  }
}

 

Suprisingly, this is just making 160 MB/s throughput. The reason is obvious when we see the disassembled code. The compiler doesn't utilize 16 xmm registers. It seems that it think there is only one xmm register. It is constantly push/pop xmm register value to the stack. This doesn't improve with -msse2 option. Or other optimization flag.

__m128i xmmTemp, xmmLeft, xmmRight, xmmEast, xmmWest, xmmCenter;
   xmmTemp = _mm_loadu_si128( (const __m128i*)(source-integerPitch-1) );
00403616  mov         edx,dword ptr [source] 
0040361C  sub         edx,dword ptr [ebp-1Ch] 
0040361F  sub         edx,1 
00403622  movdqu      xmm0,xmmword ptr [edx] 
00403626  movdqa      xmmword ptr [ebp-0D0h],xmm0 
0040362E  movdqa      xmm0,xmmword ptr [ebp-0D0h] 
00403636  movdqa      xmmword ptr [xmmTemp],xmm0 
   xmmEast = _mm_cvtepu8_epi16( xmmTemp );
0040363E  pmovzxbw    xmm0,mmword ptr [xmmTemp] 
00403647  movdqa      xmmword ptr [ebp-0F0h],xmm0 
0040364F  movdqa      xmm0,xmmword ptr [ebp-0F0h] 
00403657  movdqa      xmmword ptr [xmmEast],xmm0 
   xmmTemp = _mm_srli_si128( xmmTemp, 1 );
0040365F  movdqa      xmm0,xmmword ptr [xmmTemp] 
00403667  psrldq      xmm0,1 
0040366C  movdqa      xmmword ptr [ebp-110h],xmm0 
00403674  movdqa      xmm0,xmmword ptr [ebp-110h] 
0040367C  movdqa      xmmword ptr [xmmTemp],xmm0 

The compiler used in here is VC++ 2008 sp1. In the internet, there is saying that VC++2010 is much better than old version. But as there is intel compiler or GCC that allows inline assembler, is still a merit to pursue the intrinsics ? Wolud the compiler makes better code than hand-written code in the end ? I guess  the answer depends.  Still the fact that the VC++2008 instrincs is incompentent on SSE stands.

SSE - Pitch comparison

How fast can it be processed in the CPU side ?  We have seen that the GPU can make around 700 MB/s throughput. But what if we let it be processed in the CPU with it's all strength ?

CPU has SIMD instructions. These instructions is quite well suited with the algorithm. The pitch comparison algorithm asks the difference of each pixel. It doesn't need to know it's adjacent neighbor pixel value nor dependant to any other pixel's intermediate processing result.

Also the pitch comarision doesn't need high precision. Who cares whether the difference is 10.1 or 10.2 pixels ? The pixel value is 0 - 255 and we knows that each pixel can vary at least 1 or 2 pixels. So all we needs is at most one more point afte digit. This allows us to employ a fixed point floating point calculation.

16 bits is more than enough for the pitch comaparison. SSE 64bit provides 128 bits wides xmm registers. Each xmm registers can hold 8 words - 16 bits pixels. Each instruction can process 8 pixels at a time.

VC++ supports inline assembler on the x86 target. Here is the code.

void SSE_Pitch0( 
  const uint8_t *h_src, uint8_t *h_dst, 
  int width, int height, 
  int roiLeft, int roiTop, int roiRight, int roiBottom,
  float horPitch, float verPitch )
{
  int sizeOfDataTypeInWords = SSE::SizeOfDataType / 2;
  int integerPitch = (int)horPitch;

  uint16_t toFloor = (uint16_t)( (horPitch - (int)horPitch) * 0xFF);
  uint16_t toCeiling = 0xFF - toFloor;

  std::vector< uint16_t > toFloors( sizeOfDataTypeInWords, toFloor );
  std::vector< uint16_t > toCeilings( sizeOfDataTypeInWords, toCeiling );
  const uint16_t* toFloorsPtr = &toFloors.front();
  const uint16_t* toCeilingsPtr = &toCeilings.front();

  _asm {
    MOV     ECX, toFloorsPtr  
    MOVDQU  XMM6, [ECX]
    MOV     EDX, toCeilingsPtr
    MOVDQU  XMM7, [EDX]
  }

  for ( int y=roiTop; y<=roiBottom; y++ ) 
  {
    for ( int x=roiLeft; x<=roiRight; x+=sizeOfDataTypeInWords ) 
    {
      const uint8_t *source = h_src + width*y+x;
      uint8_t *target = h_dst + width*y+x;
      
      _asm {

        MOV      ECX, source  
        SUB      ECX, integerPitch  // east pitch pixel
        PMOVZXBW XMM1, [ECX-1]  // left
        PMOVZXBW XMM2, [ECX]    // right

        PMULLW   XMM1, XMM6  
        PMULLW   XMM2, XMM7
        PADDUSW  XMM1, XMM2  
        PSRLW    XMM1, 1

        MOV      ECX, source  
        ADD      ECX, integerPitch  // west pitch pixel

        PMOVZXBW XMM2, [ECX]    // left
        PMOVZXBW XMM3, [ECX+1]  // right

        PMULLW   XMM2, XMM7  
        PMULLW   XMM3, XMM6
        PADDUSW  XMM2, XMM3
        PSRLW    XMM2, 1

        MOV      ECX, source;  
        PMOVZXBW XMM0, [ECX]    // XMM0 = I(n) | ... | I(n+7)
        PSLLW    XMM0, 8
        MOVDQU   XMM4, XMM0

        PSUBUSW  XMM0, XMM1    // 2C - (L+R)
        PSUBUSW  XMM0, XMM2    

        PADDUSW  XMM1, XMM2    // (L+R) - 2C
        PSUBUSW  XMM1, XMM4

        PADDUSW  XMM0, XMM1

        PSRLW    XMM0, 8
        PXOR     XMM1, XMM1
        PACKUSWB XMM1, XMM0

        MOV      EDX, target
        MOVHPS   [EDX], XMM1
      }
    }
  }
}

 

Above code can be executed at around 940MB/s

CUDA Study - texture memory

CUDA provides a special memroy type so called texture memory. When image is stored in this memory, it's layout is not row-major but interleaving. This layout is supposed to give better hit and give boost on the load operation when spatially distanced pixels are used. And it provides a linear filtering in hardware which eliminate liner interpolation in the kernel.

As the memory does not have a common layout, it should be created with dedicated method - cuaMallocArray and is read-only and only accessible using special texture fetch functions.  Well in 5.0 CUDA, they say it can be written using Surface API but that's not discussed in here

The CUDA programming guide tells that there are two API for texture. One is Texture Reference API and the other is Texture Object API. The Texture Reference API is used in here. It seems that static texture variable is the only short-fall compared to Texture Object API. And Object API is introduced at 5.0 which means there is not many help if stuck.

// texture<> is only allowed as file scope static variable
// cudaReadModeElementType can't be used with cudaFilterModeLinear with uint8_t texel 
texture< uint8_t, cudaTextureType2D, cudaReadModeNormalizedFloat > gPitch2Texture;

__global__ void KernelPitch2()
{
  ...
  // add 0.5f to all texture coordinates for center offset
  float east, west, center;
  east = tex2D( gPitch2Texture, (float)( x - pc->Pitch + 0.5f), (float)(y+0.5f) ) * 255;
  west = tex2D( gPitch2Texture, (float)( x + pc->Pitch + 0.5f), (float)(y+0.5f) ) * 255;
  center = tex2D( gPitch2Texture, (float)(x+0.5f), (float)(y+0.5f) ) * 255;

  float diff = ( center * 2.0f  - ( east + west ) ) / 2.0f;
  ...
}

void CUDA_Pitch2( ... )
{
  ...
  cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc< uint8_t >();

  // texture memory should be created using cudaArray and binded with texture handle
  struct cudaArray *devSrcArray;
  cudaMallocArray( &devSrcArray, &channelDesc, width, height );
  cudaMemcpyToArray( devSrcArray, 0, 0, h_src, imgSize, cudaMemcpyHostToDevice );

  gPitch2Texture.addressMode[0] = cudaAddressModeClamp;
  gPitch2Texture.addressMode[1] = cudaAddressModeClamp;
  gPitch2Texture.filterMode = cudaFilterModeLinear;
  gPitch2Texture.normalized = false;

  cudaBindTextureToArray( gPitch2Texture, devSrcArray, channelDesc );

  const int THREAD_WIDTH = 64, THREAD_HEIGHT = 4;
  dim3 blocks( THREAD_WIDTH, THREAD_HEIGHT );
  dim3 grids( 
    (int)ceil( (double)(roiRight-roiLeft+1) / (double)blocks.x ), 
    (int)ceil( (double)(roiBottom-roiTop+1) / (double)blocks.y )
  );
 
  ...
  KernelPitch2<<< grids, blocks >>>( ... )
  ...

  cudaFreeArray( devSrcArray );
}

This makes KernelPitch2 finish in 2.117ms and makes total through put to be around 740MB/s. The last score was 630MB/s - pinned host memory and shared load version.

CUDA Study - L1 cache

Changed the thread dimensions from 16 x 16 to 64 x 4. Then the Global Memory Store Efficiency rises up from 33% to 49.9%.

Actually it rise to 49.9% at 32 x 8 and stay there no matter how large the width is. In opposite direction, if it is 8 x 32, it goes down to 20%.

This must be related with the L1 and L2 cache size. L1 cache is 16 KB and L2 256KB. L1 cache line is 128 B. When a block has 16 x 16 threads, it needs 32 ( = 16 x 2 ) cache lines - the pitch is sufficiently small like less than 64 pixels (  = 128B / 2  ). 32 cache lines corresponds to 4 KB ( = 32 cache lines x 128 B ). When maximum threads per processors is 2048, it can hold 8 blocks ( = 2048 threads / 256 threads ), and it needs 32 KB ( = 8 blocks x 4 KB ) cache to fillup all the read/write reqeust. This is exceeding 16KB limit of L1 cache. When a block is 32 x 8, it needs half and it can fit in 16KB L1 cache.

Another test is to remove the store instruction and it makes around two times fatser. What is the store efficiency ? How good can it be ? How does it being calculated ? 

The GT 640 has 0.891GB/s DDR3 with 128 Bits bus width. And it means 28.5 GB /s  ( = 0.891 x 2 x 128 / 8 ) throughput. The horizontal pitch comparison algorithm needs 5 pixels to read and 1 pixel to write.

If there is no calculation and just pure memory transaction, the maximal throughput is 4.75GB/s ( = 28.5 / 6 ). This also assumes that threads are cooperating to fully utilize the cache - 1 reads get 16 B ( = 128bits / 8 ) which will be stored in the cache and eventually consumed by 16 threads.

The current implementaion only runs around 0.4GB/s which is just 10% of maximum.

CUDA Study - shared memory

This time I am trying to use shared memory. The idea is to load source pixels collaboratively by multiple threads. If the bottle neck is the memory loading time then it will give me some boost. Though the short fall of this idea is that the last column of threads in a block should load additional pixel. Can' remove this additional load unless radically modify the whole algorithm.

Still pursuing this idea to see how it goes.

To define shared memory, it asks contant value for the size. It means blockDim.x and y can't be used for this. Though there is integer template which can conveniently declare width and height of the shared memory.

template< int W, int H >
__global__ void KernelPitch1( ... ) 
{
  ...
  __shared__ uint8_t sharedEast[ (W+1) * H ];
  __shared__ uint8_t sharedWest[ (W+1) * H ];

  int sharedIndex = threadIdx.y * (W+1) + threadIdx.x;
  sharedEast[ sharedIndex ] = d_src[ idx - pitch - 1 ];
  sharedWest[ sharedIndex ] = d_src[ idx + pitch ];

  if( threadIdx.x == (W-1) )  // one more column for each block
  {
    sharedEast[ sharedIndex + 1 ] = d_src[ idx - pitch ];
    sharedWest[ sharedIndex + 1 ] = d_src[ idx + pitch + 1 ];
  }

  __syncthreads();

  float east = 
    (float)(sharedEast[ sharedIndex   ]) * pc->ToCeiling +
    (float)(sharedEast[ sharedIndex+1 ]) * pc->ToFloor;
  float west = 
    (float)(sharedWest[ sharedIndex   ]) * pc->ToFloor + 
    (float)(sharedWest[ sharedIndex+1 ]) * pc->ToCeiling;

  float diff = ( (float)d_src[idx] * 2.0f  - ( east + west ) ) / 2.0f;
  
  if( isInBound )
  {
    d_dst[ idx ] = (unsigned char)( diff >= 0 ? diff : -diff  );
  }
}

void CUDA_Pitch1( ... )
{
  ...
  KernelPitch1< 16, 16 > <<< grids, threads >>>( ... );
  ...
}

Here is the Visual Profiler result that say slight improvement on the performance.

Though memory loading time might not be a big bottleneck in this case as the profiler points to the store operation as culprit.





One thing to notice is that the CUDA seems get confused with two level 'if' blocks. It is not a good idea to have multiple level of branches as the core has to keep all the branching but still suprised that it compiles and generate strange result. Here is the pattern that have to be avoided.

template< int W, int H >
__global__ void KernelPitch1( ... ) 
{
  ...
  
  if( isInBound )
  {
    ...

    if( threadIdx.x == (W-1) )  
    {
      // This second level if block causes trouble.
      // It seems that this block may not be executed with valid condition
      sharedEast[ sharedIndex + 1 ] = d_src[ idx - pitch ];
      sharedWest[ sharedIndex + 1 ] = d_src[ idx + pitch + 1 ];
    }

    ...
  }
}

CUDA Study - Pinned Memory

As CUDA advocates the pinned memory, only the memory allocation has been modified to see the effect. Refer below code.

struct AllocPinnedMemory {
  static uint8_t* Alloc( int width, int height ) { 
    uint8_t *p = NULL;
    cudaHostAlloc( (void**)&p, width*height, cudaHostAllocDefault );
    return p; 
  }
  static void Free( uint8_t* p ) { cudaFreeHost( p ); }
};

typedef Image < uint8_t , AllocPinnedMemory > PinnedImage;

BOOST_AUTO_TEST_CASE( TestProcess2 )
{  
  // ...
  PinnedImage src( width, height );
  PinnedImage dst( width, height );
  // ...
}

This makes the memory operation to be exactly two times faster than previous host memory which is allocated by 'new'. Refer below screenshot of the Visual Profiler. The non pinned memory copy takes 1.24 ms ( HtoD ) and 1.27 ms ( DtoH). The pinned memory takes 628 us ( both HtoD and DtoH  ). Of course there is no difference in the kernel time.

CUDA Study - Pattern Inspection Algorithm #1

There is an image processing algorithm that I want to see how fast can it be implemented. It is the first step of a common pattern inspection algorithm used in the LCD or other regular patterns.

The algorithm is simply calculate the difference between pitch distance pixels.  Here is a horizontal pitch comparison version in plain C with a fixed-point floating calculation.

uint8_t PlainC_HorizontalPitch( 
 uint8_t center,
 uint8_t eastLeft, uint8_t eastRight, uint8_t westLeft, uint8_t westRight, 
 float pitch )
{
 uint16_t toFloor = (uint16_t)( ((pitch - (int)pitch) * 0xFF)+0.5f);
 uint16_t toCeiling = 0xFF - toFloor;

 uint16_t east = toFloor * eastLeft + toCeiling * eastRight;
 uint16_t west = toFloor * westRight + toCeiling * westLeft;

 east = east >> 1;   // fx9.16
 west = west >> 1; 

 uint16_t diff1, diff2, diff;
 uint16_t c2 = center << 8;   // fx9.16   : this is 2c

 diff1 = max( 0,  c2 - east - west );
 diff2 = max( 0,  east + west - c2 );
 diff = ( diff1 + diff2 ) >> 8;   // fx8.8   : 

 return (uint8_t)diff;
} 

Now this is modified to run on the CUDA kernel like below. I really don't expect that it performs well. This version is just a very crude verions and want to see where is the bottle neck and make a plan on optimization.

__global__ 
void KernelPitch0( 
 unsigned char* d_src, unsigned char* d_dst, 
 const Area* imgSize, 
 const Region* roi,
 const PitchContext* pc )
{
 int x = roi->Left + blockIdx.x * blockDim.x + threadIdx.x;
 int y = roi->Top + blockIdx.y * blockDim.y + threadIdx.y;
 int idx = x + y * imgSize->Width;

 bool isInBound = x < roi->Right && y < roi->Bottom;

 if( isInBound )
 {
  float east = 
   (float)d_src[ idx - pc->IntegerPitch - 1 ] * pc->ToCeiling +
   (float)d_src[ idx - pc->IntegerPitch ] * pc->ToFloor;
  float west = 
   (float)d_src[ idx + pc->IntegerPitch ] * pc->ToFloor + 
   (float)d_src[ idx + pc->IntegerPitch + 1 ] * pc->ToCeiling;

  float diff = ( (float)d_src[idx] * 2.0f  - ( east + west ) ) / 2.0f;

  d_dst[ idx ] = (unsigned char)( diff >= 0 ? diff : -diff  );
 }
}

Kernel is launched as below code snippet.

void CUDA_Pitch0( 
 const unsigned char *h_src, unsigned char *h_dst, 
 int width, int height, 
 int roiLeft, int roiTop, int roiRight, int roiBottom,
 float horPitch, float verPitch )
{
 using namespace boost;

 int imgSize = width*height;

 DeviceByteImagePtr devSrc( new DeviceByteImage( width, height ) );
 DeviceByteImagePtr devDst( new DeviceByteImage( width, height ) );

 dim3 grids( width / 32, height / 32 );
 dim3 threads( 32, 32 );

 HANDLE_ERROR( cudaMemcpy( devSrc->GetPixelPtr(), h_src, imgSize, cudaMemcpyHostToDevice ) );

 shared_ptr devImageDim = CreateSmartDeviceMemoryPointer();
 shared_ptr devROI = CreateSmartDeviceMemoryPointer();
 shared_ptr devHorPitch = CreateSmartDeviceMemoryPointer();
 
 {
  Area imageDim( width, height );
  HANDLE_ERROR( cudaMemcpy( devImageDim.get(), &imageDim, sizeof(Area), cudaMemcpyHostToDevice ) );

  Region roi( roiLeft, roiTop, roiRight, roiBottom );
  HANDLE_ERROR( cudaMemcpy( devROI.get(), &roi, sizeof(Region), cudaMemcpyHostToDevice ) );

  PitchContext pc;

  pc.Pitch = horPitch;
  pc.IntegerPitch = (int)horPitch;
  pc.ToFloor = horPitch - pc.IntegerPitch;
  pc.ToCeiling = 1.0f - pc.ToFloor;
  HANDLE_ERROR( cudaMemcpy( devHorPitch.get(), &pc, sizeof(PitchContext), cudaMemcpyHostToDevice ) );
 }

 KernelPitch0<<< grids, threads >>>( 
   devSrc->GetPixelPtr(), devDst->GetPixelPtr(), 
   devImageDim.get(), devROI.get(), 
   devHorPitch.get()
  );
  if ( cudaGetLastError() != cudaSuccess ) 
  {
   cerr << "Launch Error" << endl;
   return;
  }

 HANDLE_ERROR( cudaMemcpy( h_dst, devDst->GetPixelPtr(), imgSize, cudaMemcpyDeviceToHost ) );
}

When 8192x512 ( 4MB ) image, it takes around 9291 usec which give me 428.3 MB/sec. The environment is as below.

• Geforce GT 640, CUDA 5.0
• Windows7 64 bits, Pentium G2120 @ 3.10GHz, 4GB RAM
• VS 2008

Looking into the profiling, the memory copy from host to device for source takes 1384 usec. And actual kernel takes 4081 usec. And then memory copy from devicde to host takes another 1384 usec. Refer below image. The first cudaMemcpy is the copy from host to device. The second cudaMemcpy is actually time for kernel and time for copy from device to host.