 __global__ void smem_cuda_transpose( const int m, double const * const a, double *c )
 {
 	
-/* declare a shared array */
-
-	__shared__ double smemArray[THREAD_X][THREAD_Y+1];
+/* declare a shared memory array */
+
+  __shared__ double smemArray[THREAD_X][THREAD_Y+1];
 	
-/* determine my row and column indices */
-
-	const int myRow = blockDim.x * blockIdx.x + threadIdx.x;
-	const int myCol = blockDim.y * blockIdx.y + threadIdx.y;
+/* determine my row and column indices for the error checking code */
+
+  const int myRow = blockDim.x * blockIdx.x + threadIdx.x;
+  const int myCol = blockDim.y * blockIdx.y + threadIdx.y;
 
 /* determine my row block and column block indices */
 
-	const int sourceBlockX = blockDim.x * blockIdx.x;
-	const int sourceBlockY = blockDim.y * blockIdx.y;
-
-	if( myRow < m && myCol < m )
-	{
+  const int sourceBlockX = blockDim.x * blockIdx.x;
+  const int sourceBlockY = blockDim.y * blockIdx.y;
+
+  if( myRow < m && myCol < m )
+  {
+/* read to the shared mem array */
+/* HINT: threadIdx.x should appear somewhere in the first argument to */
+/* your INDX calculation for both a[] and c[].  This will ensure proper */
+/* coalescing. */
+
+   smemArray[threadIdx.x][threadIdx.y] = 
+      a[INDX( sourceBlockX + threadIdx.x, sourceBlockY + threadIdx.y, m )];
+  } /* end if */
+
+/* synchronize */
+  __syncthreads();
 		
-		/* read to the shared mem array */
-        smemArray[threadIdx.x][threadIdx.y] = a[INDX( sourceBlockX + threadIdx.x, sourceBlockY + threadIdx.y, m )];
-		
-	} /* end if */
-		/* synchronize */
-		__syncthreads();
-		
-	if( myRow < m && myCol < m )
-	{
-		/* write the result */
-	    c[INDX( sourceBlockY + threadIdx.x, sourceBlockX + threadIdx.y, m )] = smemArray[threadIdx.y][threadIdx.x];
-	
-	} /* end if */
-	return;
+  if( myRow < m && myCol < m )
+  {
+/* write the result */
+    c[INDX( sourceBlockY + threadIdx.x, sourceBlockX + threadIdx.y, m )] = 
+           smemArray[threadIdx.y][threadIdx.x];
+  } /* end if */
+  return;
+
 } /* end naive_cuda_transpose */
 
 void host_transpose( const int m, double const * const a, double *c )

 *  limitations under the License.
 */

#include <stdio.h>
#include <math.h>


 *  naive matrix transpose goes here.
 */
 
  for( int j = 0; j < m; j++ )
  {
    for( int i = 0; i < m; i++ )
    {
      c[INDX(i,j,m)] = a[INDX(j,i,m)];
    } /* end for i */
  } /* end for j */

} /* end host_dgemm */

int main( int argc, char *argv[] )
{

  int size = SIZE;

  fprintf(stdout, "Matrix size is %d\n",size);

/* declaring pointers for array */

  double *h_a, *h_c;
  double *d_a, *d_c;
 
  size_t numbytes = (size_t) size * (size_t) size * sizeof( double );

/* allocating host memory */

  h_a = (double *) malloc( numbytes );
  if( h_a == NULL )
  {
    fprintf(stderr,"Error in host malloc h_a\n");
    return 911;
  }

  h_c = (double *) malloc( numbytes );
  if( h_c == NULL )
  {
    fprintf(stderr,"Error in host malloc h_c\n");
    return 911;
  }

/* allocating device memory */

  CUDA_CALL( cudaMalloc( (void**) &d_a, numbytes ) );
  CUDA_CALL( cudaMalloc( (void**) &d_c, numbytes ) );

/* set result matrices to zero */

  memset( h_c, 0, numbytes );
  CUDA_CALL( cudaMemset( d_c, 0, numbytes ) );

  fprintf( stdout, "Total memory required per matrix is %lf MB\n", 
     (double) numbytes / 1000000.0 );

/* initialize input matrix with random value */

  for( int i = 0; i < size * size; i++ )
  {
    h_a[i] = double( rand() ) / ( double(RAND_MAX) + 1.0 );
  } /* end for */

/* copy input matrix from host to device */

  CUDA_CALL( cudaMemcpy( d_a, h_a, numbytes, cudaMemcpyHostToDevice ) );

/* create and start timer */

  cudaEvent_t start, stop;
  CUDA_CALL( cudaEventCreate( &start ) );
  CUDA_CALL( cudaEventCreate( &stop ) );
  CUDA_CALL( cudaEventRecord( start, 0 ) );

/* call naive cpu transpose function */

  host_transpose( size, h_a, h_c );

/* stop CPU timer */

  CUDA_CALL( cudaEventRecord( stop, 0 ) );
  CUDA_CALL( cudaEventSynchronize( stop ) );
  float elapsedTime;
  CUDA_CALL( cudaEventElapsedTime( &elapsedTime, start, stop ) );

/* print CPU timing information */

  fprintf(stdout, "Total time CPU is %f sec\n", elapsedTime / 1000.0f );
  fprintf(stdout, "Performance is %f GB/s\n", 
    8.0 * 2.0 * (double) size * (double) size / 
    ( (double) elapsedTime / 1000.0 ) * 1.e-9 );

/* setup threadblock size and grid sizes */

  dim3 threads( THREAD_X, THREAD_Y, 1 );
  dim3 blocks( ( size / THREAD_X ) + 1, ( size / THREAD_Y ) + 1, 1 );

/* start timers */
  CUDA_CALL( cudaEventRecord( start, 0 ) );

/* call naive GPU transpose kernel */

  smem_cuda_transpose<<< blocks, threads >>>( size, d_a, d_c );
  CUDA_CHECK();
  CUDA_CALL( cudaDeviceSynchronize() );

/* stop the timers */

  CUDA_CALL( cudaEventRecord( stop, 0 ) );
  CUDA_CALL( cudaEventSynchronize( stop ) );
  CUDA_CALL( cudaEventElapsedTime( &elapsedTime, start, stop ) );

/* print GPU timing information */

  fprintf(stdout, "Total time GPU is %f sec\n", elapsedTime / 1000.0f );
  fprintf(stdout, "Performance is %f GB/s\n", 
    8.0 * 2.0 * (double) size * (double) size / 
    ( (double) elapsedTime / 1000.0 ) * 1.e-9 );

/* copy data from device to host */

  CUDA_CALL( cudaMemset( h_a, 0, numbytes ) );
  CUDA_CALL( cudaMemcpy( h_a, d_c, numbytes, cudaMemcpyDeviceToHost ) );

/* compare GPU to CPU for correctness */

  for( int j = 0; j < size; j++ )
  {
    for( int i = 0; i < size; i++ )
    {
      if( h_c[INDX(i,j,size)] != h_a[INDX(i,j,size)] ) 
      {
        printf("Error in element %d,%d\n", i,j );
        printf("Host %f, device %f\n",h_c[INDX(i,j,size)],
                                      h_a[INDX(i,j,size)]);
      }
    } /* end for i */
  } /* end for j */

/* free the memory */

  free( h_a );
  free( h_c );
  CUDA_CALL( cudaFree( d_a ) );
  CUDA_CALL( cudaFree( h_a ) );

  CUDA_CALL( cudaDeviceReset() );

  return 0;
}
