 
     // Create 2D Poisson problem
     cusp::coo_matrix<IndexType,ValueType,MemorySpace> A;
-//    cusp::gallery::poisson5pt(A, 10, 10);
+//    cusp::gallery::poisson5pt(A, 4, 4);
     cusp::gallery::poisson5pt(A, 50, 50);
     
     // setup linear system

        //for (int i = 0; i < levels.size(); i++)
        //    printf("level[%2d] %10d unknowns %10d nonzeros\n", i, levels[i].A.num_rows, levels[i].A.num_entries);

        //cusp::io::write_matrix_market_file(levels[0].R, "/home/nathan/Desktop/AMG/R0.mtx");
        //cusp::io::write_matrix_market_file(levels[0].A, "/home/nathan/Desktop/AMG/A0.mtx");
        //cusp::io::write_matrix_market_file(levels[0].P, "/home/nathan/Desktop/AMG/P0.mtx");
        //cusp::io::write_matrix_market_file(levels[1].A, "/home/nathan/Desktop/AMG/A1.mtx");
        //cusp::io::write_matrix_market_file(levels[1].R, "/home/nathan/Desktop/AMG/R1.mtx");
        //cusp::io::write_matrix_market_file(levels[1].P, "/home/nathan/Desktop/AMG/P1.mtx");
        //cusp::io::write_matrix_market_file(levels[2].A, "/home/nathan/Desktop/AMG/A2.mtx");

                     cusp::array1d<ValueType,cusp::device_memory>& x) const
    {
        // TODO check sizes
        const cusp::coo_matrix<IndexType,ValueType,MemorySpace> & A = levels[0].A;

        cusp::array1d<ValueType,MemorySpace> residual(A.num_rows);  // TODO eliminate temporaries
            
