+#include <cusp/csr_matrix.h>
+#include <cusp/gallery/poisson.h>
+#include <cusp/krylov/cg.h>
+#include <cusp/precond/aggregation/smoothed_aggregation.h>
 #include <cusp/relaxation/jacobi.h>
 #include <cusp/relaxation/polynomial.h>
-#include <cusp/precond/aggregation/smoothed_aggregation.h>
-#include <cusp/krylov/cg.h>
-#include <cusp/gallery/poisson.h>
-#include <cusp/csr_matrix.h>
+#include <cusp/relaxation/gauss_seidel.h>
 
 #include <iostream>
 

int main(int argc, char ** argv)
{
    typedef int                 IndexType;
    typedef float               ValueType;
    typedef cusp::device_memory MemorySpace;

    // create an empty sparse matrix structure
    cusp::hyb_matrix<IndexType, ValueType, MemorySpace> A;

    IndexType N = 256;

    // create 2D Poisson problem
    cusp::gallery::poisson5pt(A, N, N);

    {
        typedef cusp::relaxation::polynomial<ValueType,MemorySpace> Smoother;
        std::cout << "\nSolving with smoothed aggregation preconditioner and polynomial smoother" << std::endl;

        timer t0;
        cusp::precond::aggregation::smoothed_aggregation<IndexType, ValueType, MemorySpace, Smoother> M(A);
        std::cout << "constructed hierarchy in " << t0.milliseconds_elapsed() << " ms " << std::endl;

        run_amg(A,M);
    }

    // solve with smoothed aggregation algebraic multigrid preconditioner and gauss-seidel smoother
    {
        typedef cusp::relaxation::gauss_seidel<ValueType,MemorySpace> Smoother;
        std::cout << "\nSolving with smoothed aggregation preconditioner and gauss-seidel smoother" << std::endl;

        timer t0;
        cusp::precond::aggregation::smoothed_aggregation<IndexType, ValueType, MemorySpace, Smoother> M(A);
