     }
 }
 
-/* Compute the Q_n matrix used to tridiagonalize first row & colum of our block pair */
+/* Computes the Q matrix of the Householder transformation */
 __global__ void ComputeQ(float *q, float *v, int dim)
 {
     int x = threadIdx.x + blockIdx.x*blockDim.x;

#include "householder.cuh"
#include "cublas_v2.h"
#include "aux.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <cuda_runtime.h>

#define IDX2C(i,j,ld) (((j)*(ld))+(i))

using namespace std;

/* Compute each element of the Householder vector based on the thread index */
__global__ void ComputeHouseVec(float *v, float *a_col, float alpha, float r, int length)
{
    int x = threadIdx.x + blockIdx.x*blockDim.x;
    
    if(x < length && threadIdx.x != 0)
    {
        // TODO: Try with shared memory


    // TODO: Put all GPU operations into one kernel

    ComputeHouseVec<<<grid, block>>>(gpuV, gpuColumn, alpha, r, dim);
    UpdateBlockPair(block_pair, gpuV, q, dim);

    // TODO: Retreive column!

    // Copy data back to CPU and free GPU memory
    cudaMemcpy(v, gpuV, (size_t)dim*sizeof(float), cudaMemcpyDeviceToHost); CUDA_CHECK;
