     cusp::coo_matrix<int, double, cusp::host_memory> A( x.size(), g.size(), x.size()*g.n());
 
     int number = 0;
+    dg::Operator<double> forward( g.dlt().forward());
     for( unsigned i=0; i<x.size(); i++)
     {
         if (!(x[i] >= g.x0() && x[i] <= g.x1())) {

    assert( x.size() == y.size());
    cusp::coo_matrix<int, double, cusp::host_memory> A( x.size(), g.size(), x.size()*g.n()*g.n());

    int number = 0;
    for( unsigned i=0; i<x.size(); i++)
    {

        std::vector<double> px = detail::coefficients( xn, g.n()), 
                            py = detail::coefficients( yn, g.n()),
                            pz = detail::coefficients( zn, 1 );
        std::vector<double> pxF(g.n(),0), pyF(g.n(), 0);
        for( unsigned l=0; l<g.n(); l++)
            for( unsigned k=0; k<g.n(); k++)
            {
                pxF[l]+= px[k]*forward(k,l);
                pyF[l]+= py[k]*forward(k,l);
            }
        std::vector<double> pxyz( g.n()*g.n());
        for(unsigned k=0; k<g.n(); k++)
            for( unsigned j=0; j<g.n(); j++)
                pxyz[k*g.n()+j]= pz[0]*pyF[k]*pxF[j];

        //...these are the matrix coefficients with which to multiply 
        if (globalbcz == dg::DIR)
