         std::cout << " At iteration "<<i<<"\n";
         std::cout << " Error of Multigrid iterations "<<err<<"\n";
     }
+    t.toc();
+    std::cout << "Took "<<t.diff()<<"s\n";
     x = dg::evaluate( initial, grid);
+    t.tic();
     multigrid.direct_solve(multi_pol, x, b, eps);
+    t.toc();
     const double norm = dg::blas2::dot( w2d, solution);
     dg::DVec error( solution);
     dg::blas1::axpby( 1.,x,-1., solution, error);
     double err = dg::blas2::dot( w2d, error);
     err = sqrt( err/norm);
     std::cout << " Error of nested iterations "<<err<<"\n";
+    std::cout << "Took "<<t.diff()<<"s\n";
 
     int varID;
     err = nc_def_var( ncid, "x_ana", NC_DOUBLE, 2, &dim3d[1], &varID);

double sol(double x, double y)  { return sin( x)*sin(y);}
double der(double x, double y)  { return cos( x)*sin(y);}


int main()
{

    err = nc_def_var( ncid, "b_num", NC_DOUBLE, 3, dim3d, &bID);
    err = nc_def_var( ncid, "r_num", NC_DOUBLE, 3, dim3d, &rID);
    ////////////////////////////////////////////////////
    std::cout << "Type nu1 (20), nu2 (20) gamma (1) max_iter (1)\n";
    unsigned nu1, nu2, gamma;
    int max_iter = 3;
    std::cin >> nu1 >> nu2 >> gamma >> max_iter;
    x = dg::evaluate( initial, grid);
    const dg::DVec solution = dg::evaluate( sol, grid);
    t.tic();
    for( int i=0; i<max_iter; i++)
    {
        multigrid.solve(multi_pol, x, b, multi_ev, nu1, nu2, gamma, eps);
