     std::cout << " Error norm of  1d integral function "<<norm<<"\n";
     // TEST if dot throws on NaN
     dg::blas1::transform( x,x, dg::LN<double>());
-    dg::blas1::dot( x,x);
+    try{
+        dg::blas1::dot( x,x);
+    }catch ( std::exception& e)
+    {
+        std::cerr << "Error thrown as expected\n";
+        std::cerr << e.what() << std::endl;
+    }
 
     std::cout << "\nFINISHED! Continue with topology/derivatives_t.cu !\n\n";
     return 0;


    dg::Grid1d g1d( 1, 2, n, Nx);
    dg::Grid2d g2d( 1, 2, 3, 4, n, Nx, Ny);
    dg::Grid3d g3d( 1, 2, 3, 4, 5, 6, n, Nx, Ny, Nz,dg::PER,dg::PER,dg::PER);

    //test evaluation functions
    const DVec func1d = dg::construct<DVec>( dg::evaluate( exp, g1d));
    const DVec func2d = dg::construct<DVec>( dg::evaluate( function, g2d));
    const DVec func3d = dg::construct<DVec>( dg::evaluate( function, g3d));
    const DVec w1d = dg::construct<DVec>( dg::create::weights( g1d));
    const DVec w2d = dg::construct<DVec>( dg::create::weights( g2d));
    const DVec w3d = dg::construct<DVec>( dg::create::weights( g3d));
    exblas::udouble res;

    double integral = dg::blas1::dot( w1d, func1d); res.d = integral;
    std::cout << "1D integral               "<<std::setw(6)<<integral <<"\t" << res.i - 4616944842743393935  << "\n";
    double sol = (exp(2.) -exp(1));
    std::cout << "Correct integral is       "<<std::setw(6)<<sol<<std::endl;
    std::cout << "Absolute 1d error is      "<<(integral-sol)<<"\n\n";

    double integral2d = dg::blas1::dot( w2d, func2d); res.d = integral2d;
    std::cout << "2D integral               "<<std::setw(6)<<integral2d <<"\t" << res.i - 4639875759346476257<< "\n";
    double sol2d = (exp(2.)-exp(1))*(exp(4.)-exp(3));
    std::cout << "Correct integral is       "<<std::setw(6)<<sol2d<<std::endl;
    std::cout << "Absolute 2d error is      "<<(integral2d-sol2d)<<"\n\n";

    double integral3d = dg::blas1::dot( w3d, func3d); res.d = integral3d;
    std::cout << "3D integral               "<<std::setw(6)<<integral3d <<"\t" << res.i - 4675882723962622631<< "\n";
    double sol3d = sol2d*(exp(6.)-exp(5.));
    std::cout << "Correct integral is       "<<std::setw(6)<<sol3d<<std::endl;
    std::cout << "Absolute 3d error is      "<<(integral3d-sol3d)<<"\n\n";

    double norm = dg::blas2::dot( func1d, w1d, func1d); res.d = norm;
    std::cout << "Square normalized 1D norm "<<std::setw(6)<<norm<<"\t" << res.i - 4627337306989890294 <<"\n";

}
};

template<class T>
T function( T x, T y)
{
        return exp(x)*exp(y);
}
double function3d( double x, double y, double z)
{
        return exp(x)*exp(y)*exp(z);
}

int main()
{
