   datafloat *h_u   = (datafloat*) calloc(N, sizeof(datafloat));
   
   // initialize host array
-  for(int n = 0; i < N; ++i){
+  for(int n = 0;n < N; ++n){
     h_u[n] = 1;
   }
 
   // Solve discrete Laplacian
-  reduction(N, h_u);
+  sum(N, h_u);
 
   // Free the host array
   free(h_u);


  // perform reduction 10 times
  int Ntests = 10, test;
  datafloat res;

  for(test=0;test<Ntests;++test){
    reduction <<< bdim, gdim >>> (N, c_u, c_partialsum);

    // Finish reduce on host
    cudaMemcpy(h_partialsum, c_partialsum, B*sizeof(datafloat), cudaMemcpyDeviceToHost);

    datafloat psum = 0;
    for(int n=0;n<B;++n){
      psum += h_partialsum[n];
    }
    printf("sum total = %g\n", psum);

  }
  
  cudaEventRecord(endEvent, 0);
  cudaEventSynchronize(endEvent);


  datafloat *h_partialsum;

  // number of thread-blocks to partial sum u
  int GDIM = (N+BDIM-1)/BDIM;

  // allocate host array
  h_partialsum = (datafloat*) calloc(GDIM, sizeof(datafloat));

  // allocate device arrays
  cudaMalloc((void**) &c_u  , N*sizeof(datafloat));
  cudaMalloc((void**) &c_partialsum , GDIM*sizeof(datafloat));

  // copy from h_u to c_u (HOST to DEVICE)
  cudaMemcpy(c_u ,  h_u ,  N*sizeof(datafloat), cudaMemcpyHostToDevice);
