     std::vector<int> E_v;
 
     V.reserve(V_num + 1);
-    E_u.reserve(E_num);
-    E_v.reserve(E_num);
+    E.reserve(E_num);
+    E_u.reserve(E_num * 2);
+    E_v.reserve(E_num * 2);
 
     // Build V, E_u, and E_v from adjacency list representation
     int edge_index = 0;


    // Create one-dimensional blocks and grids based upon blocksize and gridsize
    dim3 dimBlock(blocksize, 1);
	dim3 dimGrid(gridsize, 1);

    // Execute CUDA kernel
    graphlets<<<dimGrid, dimBlock>>>(V_ptr, V_num, E_u_ptr, E_v_ptr, E_num, outputs_ptr);

    // Copy output data from GPU memory back into main memory
    cudaMemcpy(outputs.data(), outputs_ptr, outputs_size, cudaMemcpyDeviceToHost);

    // Free memory in GPU
	cudaFree(V_ptr);
	cudaFree(E_u_ptr);
    cudaFree(E_v_ptr);
    cudaFree(outputs_ptr);

    // Compute graphlet counts based upon edge outputs
    int g31 = 0, g32 = 0, g33 = 0, g34;

    for (int i = 0; i < E_num; i++) {
        g31 += outputs[i].tri_e;
        g32 += outputs[i].star_e;
        g33 += V_num - (outputs[i].tri_e + outputs[i].star_e + 2);
    }

    g31 /= 3;
    g32 /= 2;

    int V_num_choose_3 = (V_num * (V_num - 1) * (V_num - 2)) / (3 * 2);
    g34 = V_num_choose_3 - g31 - g32 - g33;

    std::cout << std::endl;
    std::cout << "     Graphlet Counts" << std::endl;
    std::cout << "--------------------------" << std::endl;
    std::cout << "triangle           (g31) : " << g31 << std::endl;
    std::cout << "2-star             (g32) : " << g32 << std::endl;
    std::cout << "3-node-1-edge      (g33) : " << g33 << std::endl;
    std::cout << "3-node-independent (g34) : " << g34 << std::endl;

	return EXIT_SUCCESS;
}

    // Return immediately if thread index is greater than maximum edge index
    if (i >= E_num) return;

    int ei = E[i];

    // Lookup the endpoints of the current edge
    int u = E_u[ei];
    int v = E_v[ei];

    int u_len = V[u + 1] - V[u] - 1;
    int** ptrs = new int*[u_len];
    int* lens = new int[u_len];
    int* inds = new int[u_len];

    int prev_w = -1;
    for (int edge_i = V[u], arr_i = 0; edge_i < V[u + 1]; ++edge_i && ++arr_i) {
        int w = E_v[edge_i];
        
        if (w == v) {
            arr_i--;
            continue;
        }

        ptrs[arr_i] = &E_v[V[w]];
        lens[arr_i] = V[w + 1] - V[w];
        inds[arr_i] = 0;

        while (inds[arr_i] < lens[arr_i] && ptrs[arr_i][inds[arr_i]] <= prev_w) {
            inds[arr_i]++;
        }

        prev_w = w;
    }

    int tri_e = 0, star_u = 0, star_v = 0, cliq_e = 0, cyc_e = 0;
    int iu = V[u], iv = V[v];
    while (iu < V[u + 1] || iv < V[v + 1]) {
        if (iv < V[v + 1]) {
            for (int i = 0; i < u_len; i++) {
                while (inds[i] < lens[i] && ptrs[i][inds[i]] < E_v[iv]) {
                    inds[i]++;
                }
            }
        }

        int *count_ptr = NULL;
        if (iu < V[u + 1] && (iv == V[v + 1] || E_v[iu] < E_v[iv])) {
            if (E_v[iu] == v) {
                iu++;
                continue;
            }

            star_u++;
            count_ptr = NULL;
        } else if (iu == V[u + 1] || E_v[iv] < E_v[iu]) {
            if (E_v[iv] == u) {
                iv++;
                continue;
            }

            star_v++;
            count_ptr = &cyc_e;
        } else {
            tri_e++;
            count_ptr = &cliq_e;
        }

        if (count_ptr != NULL) {
            for (int i = 0; i < u_len; i++) {
                if (inds[i] < lens[i] && ptrs[i][inds[i]] == E_v[iv]) {
                    (*count_ptr)++;
                }
            }
        }

        if (iu < V[u + 1] && (iv == V[v + 1] || E_v[iu] < E_v[iv])) {
            iu++;
        } else if (iu == V[u + 1] || E_v[iv] < E_v[iu]) {
            iv++;
        } else {
            iu++;
            iv++;
        }
    }

    outputs[i].g31 = tri_e;
    outputs[i].g32 = star_u + star_v;
    outputs[i].g33 = V_num - (tri_e + star_u + star_v + 2);
    outputs[i].g41 = cliq_e;
    outputs[i].g44 = cyc_e;

    outputs[i].T_T = (tri_e * (tri_e - 1)) / 2;
    outputs[i].Su_Sv = star_u * star_v;
    outputs[i].T_SuVSv = tri_e * (star_u + star_v);
    outputs[i].S_S = ((star_u * (star_u - 1)) / 2) + ((star_v * (star_v - 1)) / 2);

    outputs[i].T_I = tri_e * outputs[i].g33;
    outputs[i].SuVSv_I = (star_u + star_v) * outputs[i].g33;
    outputs[i].I_I = (outputs[i].g33 * (outputs[i].g33 - 1)) / 2;
    outputs[i].I_I_1 = E_num - (V[u + 1] - V[u] - 1) - (V[v + 1] - V[v] - 1) - 1;
}

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