 //     
 //     
   
-    dsNU( function, derivative); //ds(f)
+//     dsNU( function, derivative); //ds(f)
 
 //     dsNU.forward( function, derivativef); //ds(f)
 //     dsNU.backward( function, derivativeb); //ds(f)

//     std::cout << "Relative Difference is  "<< sqrt( errinvT/normf )<<"\n";
//     
    std::cout << "--------------------testing inversion dsTds" << std::endl; 
    std::cout << " # of iterations "<< invert( dsNU, functionTinv2,solutiondsTds ) << std::endl; //is dsTds
    std::cout << "Norm analytic Solution  "<<sqrt( normf)<<"\n";
    double errinvT2 =dg::blas2::dot( w3d, functionTinv2);
    std::cout << "Norm numerical Solution "<<sqrt( errinvT2)<<"\n";

//     //arithmetic average
//     dg::blas1::axpby(0.5,dsTdsb,0.5,dsTdsf,dsTdsfb);
//     dg::blas1::axpby(0.5,dsTdsbd,0.5,dsTdsfd,dsTdsfbd); 
//     ds.symv(function,dsTdsfb);
//     dg::blas1::pointwiseDot(v3d,dsTdsfb,dsTdsfb);
        //ds( function, temp);
        //dg::blas1::pointwiseDot( temp, inverseB, temp);
        //ds(temp, dsTdsfb);
        //dg::blas1::pointwiseDivide( dsTdsfb, inverseB, dsTdsfb);
//     ds.centeredT( derivative2, dsTds2); //dsT(ds(f))
//     dg::blas1::pointwiseDivide(ones,  inverseB, temp2); //B
//      ds.centeredT( ones, divbT);
//     
//     double normdsds =dg::blas2::dot(derivative2, w3d,derivative2);
//     double normds1ds =dg::blas2::dot(derivativeones, w3d,derivative2);
//     double normdivBT =dg::blas2::dot(divBT, w3d,divBT);
//      double normdivbT =dg::blas2::dot(divbT, w3d,divbT);
//      double normdivb =dg::blas2::dot(divbsol, w3d,divbsol); 
//     double normdsTf = dg::blas2::dot(derivativeT2, w3d, function2);
//     double normdsT_1 = dg::blas2::dot(derivativeT2, w3d, ones);
//     double normdsT1 = dg::blas2::dot(derivativeTones, w3d, function2);
