 
     dg::Grid1d g1d( 1, 2, n, Nx);
     dg::Grid2d g2d( 1, 2, 3, 4, n, Nx, Ny);
+    dg::RealGrid2d<float> gf2d( 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));
+    const dg::DVec func1d = dg::construct<dg::DVec>( dg::evaluate( exp, g1d));
+    const dg::DVec func2d = dg::construct<dg::DVec>( dg::evaluate( function<double>, g2d));
+    const dg::fDVec funcf2d = dg::construct<dg::fDVec>( dg::evaluate( function<float>, gf2d));
+    const dg::DVec func3d = dg::construct<dg::DVec>( dg::evaluate( function3d, g3d));
+    const dg::DVec w1d = dg::construct<dg::DVec>( dg::create::weights( g1d));
+    const dg::DVec w2d = dg::construct<dg::DVec>( dg::create::weights( g2d));
+    const dg::fDVec wf2d = dg::construct<dg::fDVec>( dg::create::weights( gf2d));
+    const dg::DVec w3d = dg::construct<dg::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";
+    std::cout << "Relative 1d error is      "<<(integral-sol)/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";
+    std::cout << "Relative 2d error is      "<<(integral2d-sol2d)/sol2d<<"\n\n";
+
+    float integralf2d = dg::blas1::dot( wf2d, func2d); res.d = integralf2d;
+    std::cout << "2D integral (float)       "<<std::setw(6)<<integralf2d <<"\t" << res.i - 4639875760323035136<< "\n";
+    float solf2d = (exp(2.)-exp(1))*(exp(4.)-exp(3));
+    std::cout << "Correct integral is       "<<std::setw(6)<<solf2d<<std::endl;
+    std::cout << "Relative 2d error (float) "<<(integralf2d-solf2d)/solf2d<<"\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";
+    std::cout << "Relative 3d error is      "<<(integral3d-sol3d)/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";

    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);

    std::cout << "\nFINISHED! Continue with topology/derivatives_t.cu !\n\n";
    return 0;

}
};

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()
{
