operator.h

Go to the documentation of this file.

/*
  This file is part of MADNESS.

  Copyright (C) 2007,2010 Oak Ridge National Laboratory

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2 of the License, or
  (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

  For more information please contact:

  Robert J. Harrison
  Oak Ridge National Laboratory
  One Bethel Valley Road
  P.O. Box 2008, MS-6367

  email: harrisonrj@ornl.gov
  tel:   865-241-3937
  fax:   865-572-0680

  $Id$
*/
#ifndef MADNESS_MRA_OPERATOR_H__INCLUDED
#define MADNESS_MRA_OPERATOR_H__INCLUDED



#include <limits.h>
#include <madness/mra/adquad.h>
#include <madness/tensor/mtxmq.h>
#include <madness/tensor/aligned.h>
#include <madness/tensor/tensor_lapack.h>
#include <madness/constants.h>

#include <madness/mra/simplecache.h>
#include <madness/mra/convolution1d.h>
#include <madness/mra/displacements.h>
#include <madness/mra/function_common_data.h>
#include <madness/mra/gfit.h>

namespace madness {


    template <typename Q, std::size_t NDIM>
    struct SeparatedConvolutionInternal {
        double norm;
        const ConvolutionData1D<Q>* ops[NDIM];
    };


    template <typename Q, std::size_t NDIM>
    struct SeparatedConvolutionData {
        std::vector< SeparatedConvolutionInternal<Q,NDIM> > muops;
        double norm;

        SeparatedConvolutionData(int rank) : muops(rank), norm(0.0) {}
        SeparatedConvolutionData(const SeparatedConvolutionData<Q,NDIM>& q) {
            muops = q.muops;
            norm = q.norm;
        }
    };



    /* this stuff is very confusing, poorly commented, and extremely poorly named!

    I think it works like this:
    We try to apply transition matrices to the compressed form of function coefficients.
    Most of the code is about caching these transition matrices. They are cached (key of the map is the displacement)
    in the SimpleCache "data", which is of type SeparatedConvolutionData, which keeps the matrices
    for all separated terms and dimensions. These SeparatedConvolutionData are constructed using
    ConvolutionND "ops", which is constructed at the construction of the SeparatedConvolution.

                        SeparatedConvolution (all terms, all dim, all displacements)

                    construction                                            storage

                                                                    SimpleCache<SeparatedConvolutionData>
                                                                        (all terms, all dim) / (all disp)
             vector<ConvolutionND>
             (1 term, all dim) / (all terms)
                                                                    vector<SeparatedConvolutionInternal>
                                                                        (1 term, all dim) / (all terms)


                                                                    vector<ConvolutionData1D>
                                                                        (1 term, 1 dim) / (all dim)

    ConvolutionND and SeparatedConvolutionInternal both point to the same data in ConvolutionData1D.
    Why we need SeparatedConvolutionInternal in the first place I have no idea. ConvolutionND has the global
    factor, and SeparatedConvolutionInternal has a norm.


    */



    template <typename Q, std::size_t NDIM>
    class SeparatedConvolution : public WorldObject< SeparatedConvolution<Q,NDIM> > {
    public:
        typedef Q opT;  
        bool doleaves;  
        bool isperiodicsum;
        bool modified_;     
        int particle_;
        bool destructive_;      

        typedef Key<NDIM> keyT;
        const static size_t opdim=NDIM;
        Timer timer_full;
        Timer timer_low_transf;
        Timer timer_low_accumulate;

        // if this is a Slater-type convolution kernel: 1-exp(-mu r12)/(2 mu)
        bool is_slaterf12;
        double mu_;

    private:


        mutable std::vector< ConvolutionND<Q,NDIM> > ops;   
        const BoundaryConditions<NDIM> bc;
        const int k;
        const FunctionCommonData<Q,NDIM>& cdata;
        const int rank;
        const std::vector<long> vk;
        const std::vector<long> v2k;
        const std::vector<Slice> s0;

        // SeparatedConvolutionData keeps data for all terms and all dimensions and 1 displacement
        mutable SimpleCache< SeparatedConvolutionData<Q,NDIM>, NDIM > data; 
        mutable SimpleCache< SeparatedConvolutionData<Q,NDIM>, 2*NDIM > mod_data; 

    public:

        bool& modified() {return modified_;}
        const bool& modified() const {return modified_;}

        int& particle() {return particle_;}
        const int& particle() const {return particle_;}

        bool& destructive() {return destructive_;}
        const bool& destructive() const {return destructive_;}

        const double& gamma() const {return mu_;}
        const double& mu() const {return mu_;}

    private:

        struct ApplyTerms {
            ApplyTerms() : r_term(false), t_term(false) {}
            bool r_term;
            bool t_term;
            bool any_terms() const {return r_term or t_term;}
        };

        struct Transformation {
            long r;             // Effective rank of transformation
            const Q* U;         // Ptr to matrix
            const Q* VT;
        };

//        /// return the right block of the upsampled operator (modified NS only)
//
//        /// unlike the operator matrices on the natural level the upsampled operator
//        /// matrices are not Toeplitz, so we need more information than just the displacement
//        ///.@param[in]  source  the source key
//        /// @param[in]  disp    the displacement
//        /// @param[in]  upop    the unfiltered operator matrix from scale n-1
//        /// @return     (k,k) patch of the upop(2k,2k) matrix
//        static Tensor<Q> operator_patch(const Translation& source, const Translation& disp, const Tensor<Q>& upop) {
//
//            // which of the 4 upsampled matrices do we need?
//            Translation sx=source%2;              // source offset
//            Translation tx=(source+disp)%2;       // target offset
//
//            Tensor<Q> rij(k,k);
//            // those two are equivalent:
//            if (sx==0 and tx==0) copy_2d_patch(rij.ptr(),             k, upop.ptr(), 2*k, k, k);
//            if (sx==1 and tx==0) copy_2d_patch(rij.ptr() + k,         k, upop.ptr(), 2*k, k, k);
//            if (sx==0 and tx==1) copy_2d_patch(rij.ptr() + 2*k*k,     k, upop.ptr(), 2*k, k, k);
//            if (sx==1 and tx==1) copy_2d_patch(rij.ptr() + 2*k*k + k, k, upop.ptr(), 2*k, k, k);
//*/
//            Slice s0(0,k-1), s1(k,2*k-1);
//            if (sx==0 and tx==0) rij=Rm(s0,s0);
//            if (sx==1 and tx==0) rij=Rm(s1,s0);
//            if (sx==0 and tx==1) rij=Rm(s0,s1);
//            if (sx==1 and tx==1) rij=Rm(s1,s1);
//
//            return rij;
//        }



        template <typename T, typename R>
        void apply_transformation(long dimk,
                                  const Transformation trans[NDIM],
                                  const Tensor<T>& f,
                                  Tensor<R>& work1,
                                  Tensor<R>& work2,
                                  Tensor<Q>& work3,
                                  const Q mufac,
                                  Tensor<R>& result) const {

            PROFILE_MEMBER_FUNC(SeparatedConvolution);
            long size = 1;
            for (std::size_t i=0; i<NDIM; ++i) size *= dimk;
            long dimi = size/dimk;

            R* restrict w1=work1.ptr();
            R* restrict w2=work2.ptr();
#ifndef HAVE_IBMBGQ
            Q* restrict w3=work3.ptr();
            const Q* U;
#endif

#ifdef HAVE_IBMBGQ
            mTxmq_padding(dimi, trans[0].r, dimk, dimk, w1, f.ptr(), trans[0].U);
#else
            U = (trans[0].r == dimk) ? trans[0].U : shrink(dimk,dimk,trans[0].r,trans[0].U,w3);
            mTxmq(dimi, trans[0].r, dimk, w1, f.ptr(), U);
#endif

            size = trans[0].r * size / dimk;
            dimi = size/dimk;
            for (std::size_t d=1; d<NDIM; ++d) {
#ifdef HAVE_IBMBGQ
                mTxmq_padding(dimi, trans[d].r, dimk, dimk, w2, w1, trans[d].U);
#else
                U = (trans[d].r == dimk) ? trans[d].U : shrink(dimk,dimk,trans[d].r,trans[d].U,w3);
                mTxmq(dimi, trans[d].r, dimk, w2, w1, U);
#endif
                size = trans[d].r * size / dimk;
                dimi = size/dimk;
                std::swap(w1,w2);
            }

            // If all blocks are full rank we can skip the transposes
            bool doit = false;
            for (std::size_t d=0; d<NDIM; ++d) doit = doit || trans[d].VT;

            if (doit) {
                for (std::size_t d=0; d<NDIM; ++d) {
                    if (trans[d].VT) {
                        dimi = size/trans[d].r;
#ifdef HAVE_IBMBGQ
                        mTxmq_padding(dimi, dimk, trans[d].r, dimk, w2, w1, trans[d].VT);
#else
                        mTxmq(dimi, dimk, trans[d].r, w2, w1, trans[d].VT);
#endif
                        size = dimk*size/trans[d].r;
                    }
                    else {
                        fast_transpose(dimk, dimi, w1, w2);
                    }
                    std::swap(w1,w2);
                }
            }
            // Assuming here that result is contiguous and aligned
            aligned_axpy(size, result.ptr(), w1, mufac);
        }


        template <typename T, typename R>
        void apply_transformation3(const Tensor<T> trans2[NDIM],
                                  const Tensor<T>& f,
                                  const Q mufac,
                                  Tensor<R>& result) const {

            PROFILE_MEMBER_FUNC(SeparatedConvolution);

            Tensor<R> result2=general_transform(f,trans2);
            result2.scale(mufac);
            result+=result2;

        }



        template <typename T, typename R>
        void apply_transformation2(Level n, long dimk,  double tol,
                                  const Tensor<T> trans2[NDIM],
                                  const GenTensor<T>& f,
                                  GenTensor<R>& work1,
                                  GenTensor<R>& work2,
                                  GenTensor<Q>& work3,
                                  const Q mufac,
                                  GenTensor<R>& result) const {

            PROFILE_MEMBER_FUNC(SeparatedConvolution);

#if 1
            result=general_transform(f,trans2);
            result.scale(mufac);

#else

            long size = 1;
            for (std::size_t i=0; i<NDIM; ++i) size *= dimk;
            long dimi = size/dimk;

            R* restrict w1=work1.ptr();
            R* restrict w2=work2.ptr();
            Q* restrict w3=work3.ptr();

            const Q* U;

            U = (trans[0].r == dimk) ? trans[0].U : shrink(dimk,dimk,trans[0].r,trans[0].U,w3);
            mTxmq(dimi, trans[0].r, dimk, w1, f.ptr(), U);
            size = trans[0].r * size / dimk;
            dimi = size/dimk;
            for (std::size_t d=1; d<NDIM; ++d) {
                U = (trans[d].r == dimk) ? trans[d].U : shrink(dimk,dimk,trans[d].r,trans[d].U,w3);
                mTxmq(dimi, trans[d].r, dimk, w2, w1, U);
                size = trans[d].r * size / dimk;
                dimi = size/dimk;
                std::swap(w1,w2);
            }

            // If all blocks are full rank we can skip the transposes
            bool doit = false;
            for (std::size_t d=0; d<NDIM; ++d) doit = doit || trans[d].VT;

            if (doit) {
                for (std::size_t d=0; d<NDIM; ++d) {
                    if (trans[d].VT) {
                        dimi = size/trans[d].r;
                        mTxmq(dimi, dimk, trans[d].r, w2, w1, trans[d].VT);
                        size = dimk*size/trans[d].r;
                    }
                    else {
                        fast_transpose(dimk, dimi, w1, w2);
                    }
                    std::swap(w1,w2);
                }
            }
            // Assuming here that result is contiguous and aligned
            aligned_axpy(size, result.ptr(), w1, mufac);
            //    long one = 1;
            //daxpy_(&size, &mufac, w1, &one, result.ptr(), &one);
#endif
        }


        template <typename T>
        void muopxv_fast(ApplyTerms at,
                         const ConvolutionData1D<Q>* const ops_1d[NDIM],
                         const Tensor<T>& f, const Tensor<T>& f0,
                         Tensor<TENSOR_RESULT_TYPE(T,Q)>& result,
                         Tensor<TENSOR_RESULT_TYPE(T,Q)>& result0,
                         double tol,
                         const Q mufac,
                         Tensor<TENSOR_RESULT_TYPE(T,Q)>& work1,
                         Tensor<TENSOR_RESULT_TYPE(T,Q)>& work2,
                         Tensor<Q>& work5) const {

            PROFILE_MEMBER_FUNC(SeparatedConvolution);
            Transformation trans[NDIM];
            Tensor<T> trans2[NDIM];

            double Rnorm = 1.0;
            for (std::size_t d=0; d<NDIM; ++d) Rnorm *= ops_1d[d]->Rnorm;

            if (at.r_term and (Rnorm > 1.e-20)) {

                tol = tol/(Rnorm*NDIM);  // Errors are relative within here

                // Determine rank of SVD to use or if to use the full matrix
                long twok = 2*k;
                if (modified()) twok=k;

                long break_even;
                if (NDIM==1) break_even = long(0.5*twok);
                else if (NDIM==2) break_even = long(0.6*twok);
                else if (NDIM==3) break_even=long(0.65*twok);
                else break_even=long(0.7*twok);
                for (std::size_t d=0; d<NDIM; ++d) {
                    long r;
                    for (r=0; r<twok; ++r) {
                        if (ops_1d[d]->Rs[r] < tol) break;
                    }
                    if (r >= break_even) {
                        trans[d].r = twok;
                        trans[d].U = ops_1d[d]->R.ptr();
                        trans[d].VT = 0;
                    }
                    else {
                        r += (r&1L);
                        trans[d].r = std::max(2L,r);
                        trans[d].U = ops_1d[d]->RU.ptr();
                        trans[d].VT = ops_1d[d]->RVT.ptr();
                    }
                    trans2[d]=ops_1d[d]->R;
                }
                apply_transformation(twok, trans, f, work1, work2, work5, mufac, result);
    //            apply_transformation2(n, twok, tol, trans2, f, work1, work2, work5, mufac, result);
//                apply_transformation3(trans2, f, mufac, result);
            }

            double Tnorm = 1.0;
            for (std::size_t d=0; d<NDIM; ++d) Tnorm *= ops_1d[d]->Tnorm;

            if (at.t_term and (Tnorm>0.0)) {
                tol = tol/(Tnorm*NDIM);  // Errors are relative within here

                long break_even;
                if (NDIM==1) break_even = long(0.5*k);
                else if (NDIM==2) break_even = long(0.6*k);
                else if (NDIM==3) break_even=long(0.65*k);
                else break_even=long(0.7*k);
                for (std::size_t d=0; d<NDIM; ++d) {
                    long r;
                    for (r=0; r<k; ++r) {
                        if (ops_1d[d]->Ts[r] < tol) break;
                    }
                    if (r >= break_even) {
                        trans[d].r = k;
                        trans[d].U = ops_1d[d]->T.ptr();
                        trans[d].VT = 0;
                    }
                    else {
                        r += (r&1L);
                        trans[d].r = std::max(2L,r);
                        trans[d].U = ops_1d[d]->TU.ptr();
                        trans[d].VT = ops_1d[d]->TVT.ptr();
                    }
                    trans2[d]=ops_1d[d]->T;
                }
                apply_transformation(k, trans, f0, work1, work2, work5, -mufac, result0);
//                apply_transformation2(n, k, tol, trans2, f0, work1, work2, work5, -mufac, result0);
//                apply_transformation3(trans2, f0, -mufac, result0);
            }
        }


        template <typename T>
        void muopxv_fast2(Level n,
                         const ConvolutionData1D<Q>* const ops_1d[NDIM],
                         const GenTensor<T>& f, const GenTensor<T>& f0,
                         GenTensor<TENSOR_RESULT_TYPE(T,Q)>& result,
                         GenTensor<TENSOR_RESULT_TYPE(T,Q)>& result0,
                         double tol,
                         const Q mufac,
                         GenTensor<TENSOR_RESULT_TYPE(T,Q)>& work1,
                         GenTensor<TENSOR_RESULT_TYPE(T,Q)>& work2,
                         GenTensor<Q>& work5) const {

            PROFILE_MEMBER_FUNC(SeparatedConvolution);
//            Transformation trans[NDIM];
            Tensor<T> trans2[NDIM];
//            MADNESS_EXCEPTION("no muopxv_fast2",1);

            double Rnorm = 1.0;
            for (std::size_t d=0; d<NDIM; ++d) Rnorm *= ops_1d[d]->Rnorm;
            if (Rnorm == 0.0) return;

            if (Rnorm > 1.e-20) {

                                tol = tol/(Rnorm*NDIM);  // Errors are relative within here

                                // Determine rank of SVD to use or if to use the full matrix
                                long twok = 2*k;
                                if (modified()) twok=k;
//                              long break_even;
//                              if (NDIM==1) break_even = long(0.5*twok);
//                              else if (NDIM==2) break_even = long(0.6*twok);
//                              else if (NDIM==3) break_even=long(0.65*twok);
//                              else break_even=long(0.7*twok);
                                for (std::size_t d=0; d<NDIM; ++d) {
                                        long r;
                                        for (r=0; r<twok; ++r) {
                                                if (ops_1d[d]->Rs[r] < tol) break;
                                        }
//                                      if (r >= break_even) {
//                                              trans[d].r = twok;
//                                              trans[d].U = ops_1d[d]->R.ptr();
//                                              trans[d].VT = 0;
//                                      }
//                                      else {
//                                              r += (r&1L);
//                                              trans[d].r = std::max(2L,r);
//                                              trans[d].U = ops_1d[d]->RU.ptr();
//                                              trans[d].VT = ops_1d[d]->RVT.ptr();
//                                      }
                                        trans2[d]=ops_1d[d]->R;
                                }
                                apply_transformation2(n, twok, tol, trans2, f, work1, work2, work5, mufac, result);
            }

            double Tnorm = 1.0;
            for (std::size_t d=0; d<NDIM; ++d) Tnorm *= ops_1d[d]->Tnorm;

            if (n > 0 and (Tnorm>1.e-20)) {
//                              long break_even;
//
//                if (NDIM==1) break_even = long(0.5*k);
//                else if (NDIM==2) break_even = long(0.6*k);
//                else if (NDIM==3) break_even=long(0.65*k);
//                else break_even=long(0.7*k);
                for (std::size_t d=0; d<NDIM; ++d) {
                    long r;
                    for (r=0; r<k; ++r) {
                        if (ops_1d[d]->Ts[r] < tol) break;
                    }
//                    if (r >= break_even) {
//                        trans[d].r = k;
//                        trans[d].U = ops_1d[d]->T.ptr();
//                        trans[d].VT = 0;
//                    }
//                    else {
//                        r += (r&1L);
//                        trans[d].r = std::max(2L,r);
//                        trans[d].U = ops_1d[d]->TU.ptr();
//                        trans[d].VT = ops_1d[d]->TVT.ptr();
//                    }
                    trans2[d]=ops_1d[d]->T;
                }
                apply_transformation2(n, k, tol, trans2, f0, work1, work2, work5, -mufac, result0);
            }
        }


        double munorm2(Level n, const ConvolutionData1D<Q>* ops[]) const {
            if (modified()) return munorm2_modified(n,ops);
            return munorm2_ns(n,ops);
        }

        double munorm2_ns(Level n, const ConvolutionData1D<Q>* ops[]) const {
            PROFILE_MEMBER_FUNC(SeparatedConvolution);
            
            double prodR=1.0, prodT=1.0;
            for (std::size_t d=0; d<NDIM; ++d) {
                prodR *= ops[d]->Rnormf;
                prodT *= ops[d]->Tnormf;

            }
            if (n) prodR = sqrt(std::max(prodR*prodR - prodT*prodT,0.0));

            // this kicks in if the line above has no numerically significant digits.
            if (prodR < 1e-8*prodT) {
                double prod=1.0, sum=0.0;
                for (std::size_t d=0; d<NDIM; ++d) {
                    double a = ops[d]->NSnormf;
                    double b = ops[d]->Tnormf;
                    double aa = std::min(a,b);
                    double bb = std::max(a,b);
                    prod *= bb;
                    if (bb > 0.0) sum +=(aa/bb);
                }
                if (n) prod *= sum;
                prodR = prod;
            }

            return prodR;
        }


        double munorm2_modified(Level n, const ConvolutionData1D<Q>* ops_1d[]) const {
            PROFILE_MEMBER_FUNC(SeparatedConvolution);

            // follows Eq. (21) ff of Beylkin 2008 (Beylkin Appl. Comput. Harmon. Anal. 24, pp 354)

            // we have all combinations of difference, upsampled, F terms (d, u, f),
            // with the constraint that d is in each term exactly once. In the mixed terms (udf)
            // we just get all possible combinations, in the pure terms (dff, duu) we have
            // to multiply each term (dff, fdf, ffd) with (NDIM-1)!, to get the right number.

            double dff = 0.0;
            double duu = 0.0;
            double udf = 0.0;

            // loop over d shifting over the dimensions dxx, xdx, xxd,
            for (size_t d=0; d<NDIM; ++d) {
                double dff_tmp = ops_1d[d]->N_diff;
                double duu_tmp = ops_1d[d]->N_diff;
                double udf_tmp = ops_1d[d]->N_diff;


                for (size_t dd=0; dd<NDIM; ++dd) {
                    if (dd!=d) {
                        dff_tmp *= ops_1d[dd]->N_F;
                        duu_tmp *= ops_1d[dd]->N_up;

                        udf_tmp *= ops_1d[dd]->N_F;
                        for (size_t ddd=0; ddd<NDIM; ++ddd) {
                            if (ddd!=dd) udf += udf_tmp * ops_1d[ddd]->N_up;
                        }
                    }
                }

                dff+=dff_tmp;
                duu+=duu_tmp;
            }

            // finalize with the factorial
            double factorial=1.0;
            for (int i=1; i<static_cast<int>(NDIM)-1; ++i) factorial*=double(i);
            dff*=factorial;
            duu*=factorial;

            // Eq. (23) of Beylkin 2008, for one separated term WITHOUT the factor
            double norm=(dff + udf + duu) /(factorial * double(NDIM));

//            // double check
//            if (NDIM==3) {
//                Tensor<Q> R_full=outer(ops_1d[0]->R,outer(ops_1d[1]->R,ops_1d[2]->R));
//                Tensor<Q> T_full=outer(ops_1d[0]->T,outer(ops_1d[1]->T,ops_1d[2]->T));
//                double n2=(R_full-T_full).normf();
//                norm=n2;
//            }

            return norm;
            
        }



        const SeparatedConvolutionInternal<Q,NDIM> getmuop(int mu, Level n, const Key<NDIM>& disp) const {
            PROFILE_MEMBER_FUNC(SeparatedConvolution);
            SeparatedConvolutionInternal<Q,NDIM> op;
            for (std::size_t d=0; d<NDIM; ++d) {
                op.ops[d] = ops[mu].getop(d)->nonstandard(n, disp.translation()[d]);
            }
            op.norm = munorm2(n, op.ops)*std::abs(ops[mu].getfac());

//             double newnorm = munorm2(n, op.ops);
//             // This rescaling empirically based upon BSH separated expansion
//             // ... needs more testing.  OK also for TDSE.
//             // All is good except for some 000 blocks which are up to sqrt(k^d) off.
//             for (int d=0; d<NDIM; ++d)  {
//                 if (disp[d] == 0) newnorm *= 0.5;
//                 else if (std::abs(disp[d]) == 1) newnorm *= 0.8;
//             }
//            double oldnorm = munorm(n, op.ops);
//             if (oldnorm > 1e-13 && (newnorm < 0.5*oldnorm || newnorm > 2.0*oldnorm) )
//                 print("munorm", n, disp, mu, newnorm, oldnorm, newnorm/oldnorm);

            return op;
        }




        const SeparatedConvolutionInternal<Q,NDIM>
        getmuop_modified(int mu, Level n, const Key<NDIM>& disp, const Key<NDIM>& source) const {
            PROFILE_MEMBER_FUNC(SeparatedConvolution);


            // SeparatedConvolutionInternal keeps data for 1 term and all dimensions
            SeparatedConvolutionInternal<Q,NDIM> op;

            // in the modified NS form we need not only the displacement, but also the source Translation
            // for correctly constructing the operator, b/c the operator is not Toeplitz

            // op.ops is of type ConvolutionData1D (1 term, 1 dim, 1 disp)
            // ops[mu] is of type ConvolutionND (1 term, all dim, 1 disp)
            for (std::size_t d=0; d<NDIM; ++d) {
                Translation sx=source.translation()[d];                          // source translation
                Translation tx=source.translation()[d]+disp.translation()[d];    // target translation

                Key<2> op_key(n,Vector<Translation,2>(vec(sx,tx)));
                op.ops[d] = ops[mu].getop(d)->mod_nonstandard(op_key);
            }

            // works for both modified and not modified NS form
            op.norm = munorm2(n, op.ops)*std::abs(ops[mu].getfac());
//            op.norm=1.0;
            return op;
        }

        const SeparatedConvolutionData<Q,NDIM>* getop(Level n, const Key<NDIM>& d, const Key<NDIM>& source) const {

            // in the NS form the operator depends only on the displacement
            if (not modified()) return getop_ns(n,d);
            return getop_modified(n, d, source);
        }



        const SeparatedConvolutionData<Q,NDIM>* getop_ns(Level n, const Key<NDIM>& d) const {
            PROFILE_MEMBER_FUNC(SeparatedConvolution);
            const SeparatedConvolutionData<Q,NDIM>* p = data.getptr(n,d);
            if (p) return p;

            // get the data for each term
            SeparatedConvolutionData<Q,NDIM> op(rank);
            for (int mu=0; mu<rank; ++mu) {
                // op.muops is of type SeparatedConvolutionInternal (1 term, all dim, 1 disp)
                // getmuop uses ConvolutionND
                op.muops[mu] = getmuop(mu, n, d);
            }

            double norm = 0.0;
            for (int mu=0; mu<rank; ++mu) {
                const double munorm = op.muops[mu].norm;
                norm += munorm*munorm;
            }
            //print("getop", n, d, norm);
            op.norm = sqrt(norm);
            data.set(n, d, op);
            return data.getptr(n,d);
        }




        const SeparatedConvolutionData<Q,NDIM>* getop_modified(Level n, const Key<NDIM>& disp, const Key<NDIM>& source) const {
            PROFILE_MEMBER_FUNC(SeparatedConvolution);

            // in the modified NS form the upsampled part of the operator depends on the modulus of the source
            Vector<Translation,NDIM> t=source.translation();
            for (size_t i=0; i<NDIM; ++i) t[i]=t[i]%2;
            Key<2*NDIM> key=disp.merge_with(Key<NDIM>(source.level(),t));

            const SeparatedConvolutionData<Q,NDIM>* p = mod_data.getptr(n,key);
            if (p) return p;

            // get the data for each term
            // op.muops is of type SeparatedConvolutionInternal (1 term, all dim, 1 disp)
            // getmuop uses ConvolutionND
            SeparatedConvolutionData<Q,NDIM> op(rank);
            for (int mu=0; mu<rank; ++mu) op.muops[mu] = getmuop_modified(mu, n, disp, source);

            double norm = 0.0;
            for (int mu=0; mu<rank; ++mu) {
                const double munorm = op.muops[mu].norm;
                norm += munorm*munorm;
            }

            op.norm = sqrt(norm);
            mod_data.set(n, key, op);
            return mod_data.getptr(n,key);
        }


        void check_cubic() {
            // !!! NB ... cell volume obtained from global defaults
            const Tensor<double>& cell_width = FunctionDefaults<NDIM>::get_cell_width();
            // Check that the cell is cubic since currently is assumed
            for (std::size_t d=1; d<NDIM; ++d) {
                MADNESS_ASSERT(fabs(cell_width(d)-cell_width(0L)) < 1e-14*cell_width(0L));
            }
        }



        template<typename T, size_t FDIM>
        GenTensor<T> partial_upsample(const Key<FDIM>& key, const GenTensor<T>& coeff, const int particle) const {

            if (coeff.rank()==0) return GenTensor<T>();
            MADNESS_ASSERT(coeff.dim(0)==k);
            if (NDIM==coeff.ndim()) {
                MADNESS_ASSERT(particle==1);    // other particle, leave this particle unchanged
                return coeff;
            }

            MADNESS_ASSERT(coeff.ndim()==FDIM);
            MADNESS_ASSERT(particle==0 or (2*NDIM==FDIM));

            // the twoscale coefficients: for upsampling use h0/h1; see Alpert Eq (3.35a/b)
            // handle the spectator dimensions with the identity matrix
            const Tensor<T> h[2] = {cdata.h0, cdata.h1};
            Tensor<T> identity(k,k);
            for (int i=0; i<k; ++i) identity(i,i)=1.0;
            Tensor<T> matrices[2*NDIM];

            // get the appropriate twoscale coefficients for each dimension
            if (particle==0) {
                for (size_t ii=0; ii<NDIM; ++ii) matrices[ii]=h[key.translation()[ii]%2];
                for (size_t ii=0; ii<NDIM; ++ii) matrices[ii+NDIM]=identity;
            } else if (particle==1) {
                for (size_t ii=0; ii<NDIM; ++ii) matrices[ii]=identity;
                for (size_t ii=0; ii<NDIM; ++ii) matrices[ii+NDIM]=h[key.translation()[ii+NDIM]%2];
            } else {
                MADNESS_EXCEPTION("unknown particle",1);
            }

            // transform and accumulate on the result
            const GenTensor<T> result=general_transform(coeff,matrices);
            return result;
        }



        template<typename T, size_t FDIM>
        GenTensor<T> upsample(const Key<FDIM>& key, const GenTensor<T>& coeff) const {

            // the twoscale coefficients: for upsampling use h0/h1; see Alpert Eq (3.35a/b)
            // note there are no difference coefficients; if you want that use unfilter
            const Tensor<T> h[2] = {cdata.h0, cdata.h1};
            Tensor<T> matrices[FDIM];

            // get the appropriate twoscale coefficients for each dimension
            for (size_t ii=0; ii<FDIM; ++ii) matrices[ii]=h[key.translation()[ii]%2];

            // transform and accumulate on the result
            const GenTensor<T> result=general_transform(coeff,matrices);
            return result;
        }


    public:

        // For separated convolutions with same operator in each direction (isotropic)
        SeparatedConvolution(World& world,
                             std::vector< std::shared_ptr< Convolution1D<Q> > >& argops,
                             const BoundaryConditions<NDIM>& bc = FunctionDefaults<NDIM>::get_bc(),
                             long k = FunctionDefaults<NDIM>::get_k(),
                             bool doleaves = false)
                : WorldObject< SeparatedConvolution<Q,NDIM> >(world)
                , doleaves(doleaves)
                , isperiodicsum(bc(0,0)==BC_PERIODIC)
                , modified_(false)
                , particle_(1)
                , destructive_(false)
                , is_slaterf12(false)
                , mu_(0.0)
                , bc(bc)
                , k(k)
                , cdata(FunctionCommonData<Q,NDIM>::get(k))
                , rank(argops.size())
                , vk(NDIM,k)
                , v2k(NDIM,2*k)
                , s0(std::max<std::size_t>(2,NDIM),Slice(0,k-1))
        {
            // Presently we must have periodic or non-periodic in all dimensions.
            for (std::size_t d=1; d<NDIM; ++d) {
                MADNESS_ASSERT(bc(d,0)==bc(0,0));
            }
            check_cubic();

            for (unsigned int mu=0; mu < argops.size(); ++mu) {
              this->ops.push_back(ConvolutionND<Q,NDIM>(argops[mu]));
            }

            this->process_pending();
        }

        // For general convolutions
        SeparatedConvolution(World& world,
                             std::vector< ConvolutionND<Q,NDIM> >& argops,
                             const BoundaryConditions<NDIM>& bc = FunctionDefaults<NDIM>::get_bc(),
                             long k = FunctionDefaults<NDIM>::get_k(),
                             bool doleaves = false)
                : WorldObject< SeparatedConvolution<Q,NDIM> >(world)
                , doleaves(doleaves)
                , isperiodicsum(bc(0,0)==BC_PERIODIC)
                , modified_(false)
                , particle_(1)
                , destructive_(false)
                , is_slaterf12(false)
                , mu_(0.0)
                , ops(argops)
                , bc(bc)
                , k(k)
                , cdata(FunctionCommonData<Q,NDIM>::get(k))
                , rank(argops.size())
                , vk(NDIM,k)
                , v2k(NDIM,2*k)
                , s0(std::max<std::size_t>(2,NDIM),Slice(0,k-1))
        {
            // Presently we must have periodic or non-periodic in all dimensions.
            for (std::size_t d=1; d<NDIM; ++d) {
                MADNESS_ASSERT(bc(d,0)==bc(0,0));
            }
            this->process_pending();
        }

        SeparatedConvolution(World& world,
                             const Tensor<Q>& coeff, const Tensor<double>& expnt,
                             const BoundaryConditions<NDIM>& bc = FunctionDefaults<NDIM>::get_bc(),
                             int k=FunctionDefaults<NDIM>::get_k(),
                             bool doleaves = false,
                             double mu=0.0)
                : WorldObject< SeparatedConvolution<Q,NDIM> >(world)
                , doleaves(doleaves)
                , isperiodicsum(bc(0,0)==BC_PERIODIC)
                , modified_(false)
                , particle_(1)
                , destructive_(false)
                , is_slaterf12(mu>0.0)
                , mu_(mu)
                , ops(coeff.dim(0))
                , bc(bc)
                , k(k)
                , cdata(FunctionCommonData<Q,NDIM>::get(k))
                , rank(coeff.dim(0))
                , vk(NDIM,k)
                , v2k(NDIM,2*k)
                , s0(std::max<std::size_t>(2,NDIM),Slice(0,k-1))
        {
            // Presently we must have periodic or non-periodic in all dimensions.
            for (std::size_t d=1; d<NDIM; ++d) {
                MADNESS_ASSERT(bc(d,0)==bc(0,0));
            }

            const Tensor<double>& width = FunctionDefaults<NDIM>::get_cell_width();
            const double pi = constants::pi;

            for (int mu=0; mu<rank; ++mu) {
                Q c = std::pow(sqrt(expnt(mu)/pi),static_cast<int>(NDIM)); // Normalization coeff

                // We cache the normalized operator so the factor is the value we must multiply
                // by to recover the coeff we want.
                ops[mu].setfac(coeff(mu)/c);

                for (std::size_t d=0; d<NDIM; ++d) {
                  ops[mu].setop(d,GaussianConvolution1DCache<Q>::get(k, expnt(mu)*width[d]*width[d], 0, isperiodicsum));
                }
            }
        }

        SeparatedConvolution(World& world,
                             Vector<double,NDIM> args,
                             const Tensor<Q>& coeff, const Tensor<double>& expnt,
                             const BoundaryConditions<NDIM>& bc = FunctionDefaults<NDIM>::get_bc(),
                             int k=FunctionDefaults<NDIM>::get_k(),
                             bool doleaves=false)
                : WorldObject< SeparatedConvolution<Q,NDIM> >(world)
                , doleaves(doleaves)
                , isperiodicsum(bc(0,0)==BC_PERIODIC)
                , modified_(false)
                , particle_(1)
                , destructive_(false)
                , is_slaterf12(false)
                , mu_(0.0)
                , ops(coeff.dim(0))
                , bc(bc)
                , k(k)
                , cdata(FunctionCommonData<Q,NDIM>::get(k))
                , rank(coeff.dim(0))
                , vk(NDIM,k)
                , v2k(NDIM,2*k)
                , s0(std::max<std::size_t>(2,NDIM),Slice(0,k-1))
        {
            // Presently we must have periodic or non-periodic in all dimensions.
            for (std::size_t d=1; d<NDIM; ++d) {
                MADNESS_ASSERT(bc(d,0)==bc(0,0));
            }

            const Tensor<double>& width = FunctionDefaults<NDIM>::get_cell_width();

            for (int mu=0; mu<rank; ++mu) {
                for (std::size_t d=0; d<NDIM; ++d) {
                  Q c = pow(coeff[mu],1.0/NDIM);
                  std::shared_ptr<GaussianConvolution1D<Q> >
                      gcptr(new GaussianConvolution1D<Q>(k, c*width[d], expnt(mu)*width[d]*width[d],
                              0, isperiodicsum, args[d]));
                  ops[mu].setop(d,gcptr);
                }
            }
        }

        virtual ~SeparatedConvolution() { }

        void print_timer() const {
                if (this->world.rank()==0) {
                timer_full.print("op full tensor       ");
                timer_low_transf.print("op low rank transform");
                timer_low_accumulate.print("op low rank addition ");
                }
        }

        void reset_timer() const {
                if (this->world.rank()==0) {
                timer_full.reset();
                timer_low_transf.reset();
                timer_low_accumulate.reset();
                }
        }

        const BoundaryConditions<NDIM>& get_bc() const {return bc;}

        const std::vector< Key<NDIM> >& get_disp(Level n) const {
            return Displacements<NDIM>().get_disp(n, isperiodicsum);
        }

        double norm(Level n, const Key<NDIM>& d, const Key<NDIM>& source_key) const {
            // SeparatedConvolutionData keeps data for all terms and all dimensions and 1 displacement
//            return 1.0;
            return getop(n, d, source_key)->norm;
        }


        template<size_t FDIM>
        typename disable_if_c<FDIM==NDIM, Key<NDIM> >::type
        get_source_key(const Key<FDIM> key) const {
            Key<NDIM> source;
            Key<FDIM-NDIM> dummykey;
            if (particle()==1) key.break_apart(source,dummykey);
            if (particle()==2) key.break_apart(dummykey,source);
                return source;
        }


        template<size_t FDIM>
        typename enable_if_c<FDIM==NDIM, Key<NDIM> >::type
        get_source_key(const Key<FDIM> key) const {
                return key;
        }


        template <typename T, size_t FDIM>
        Function<TENSOR_RESULT_TYPE(T,Q),FDIM> operator()(const Function<T,FDIM>& f) const {
            return madness::apply(*this, f);
        }


        template <typename T, size_t LDIM>
        Function<TENSOR_RESULT_TYPE(T,Q),LDIM+LDIM>
        operator()(const Function<T,LDIM>& f1, const Function<Q,LDIM>& f2) const {
                MADNESS_ASSERT(not is_slaterf12);
            return madness::apply(*this, f1, f2);
        }



        template <typename T>
        Tensor<TENSOR_RESULT_TYPE(T,Q)> apply(const Key<NDIM>& source,
                                              const Key<NDIM>& shift,
                                              const Tensor<T>& coeff,
                                              double tol) const {
            PROFILE_MEMBER_FUNC(SeparatedConvolution);
            MADNESS_ASSERT(coeff.ndim()==NDIM);

            double cpu0=cpu_time();

            typedef TENSOR_RESULT_TYPE(T,Q) resultT;
            const Tensor<T>* input = &coeff;
            Tensor<T> dummy;

            if (not modified()) {
                if (coeff.dim(0) == k) {
                    // This processes leaf nodes with only scaling
                    // coefficients ... FuncImpl::apply by default does not
                    // apply the operator to these since for smoothing operators
                    // it is not necessary.  It is necessary for operators such
                    // as differentiation and time evolution and will also occur
                    // if the application of the operator widens the tree.
                    dummy = Tensor<T>(v2k);
                    dummy(s0) = coeff;
                    input = &dummy;
                }
                else {
                    MADNESS_ASSERT(coeff.dim(0)==2*k);
                }
            }

            tol = tol/rank; // Error is per separated term
            ApplyTerms at;
            at.r_term=true;
            at.t_term=(source.level()>0);

            const SeparatedConvolutionData<Q,NDIM>* op = getop(source.level(), shift, source);

            //print("sepop",source,shift,op->norm,tol);

            Tensor<resultT> r(v2k), r0(vk);
            Tensor<resultT> work1(v2k,false), work2(v2k,false);
            Tensor<Q> work5(2*k,2*k);

            if (modified()) {
                   r=Tensor<resultT>(vk);
                   work1=Tensor<resultT>(vk,false);
                   work2=Tensor<resultT>(vk,false);
                   work5=Tensor<Q>(k,k);
            }

            const Tensor<T> f0 = copy(coeff(s0));
            for (int mu=0; mu<rank; ++mu) {
                // SeparatedConvolutionInternal keeps data for 1 term and all dimensions and 1 displacement
                const SeparatedConvolutionInternal<Q,NDIM>& muop =  op->muops[mu];
                if (muop.norm > tol) {
                    // ops is of ConvolutionND, returns data for 1 term and all dimensions
                    Q fac = ops[mu].getfac();
                    muopxv_fast(at, muop.ops, *input, f0, r, r0, tol/std::abs(fac), fac,
                                work1, work2, work5);
                }
            }

            r(s0).gaxpy(1.0,r0,1.0);
            double cpu1=cpu_time();
            timer_full.accumulate(cpu1-cpu0);

            return r;
        }



        template<typename T>
        GenTensor<TENSOR_RESULT_TYPE(T,Q)> apply2_lowdim(const Key<NDIM>& source,
                const Key<NDIM>& shift, const GenTensor<T>& coeff, double tol, double tol2) const {

            typedef TENSOR_RESULT_TYPE(T,Q) resultT;

            // some checks
            MADNESS_ASSERT(coeff.tensor_type()==TT_2D);           // for now
            MADNESS_ASSERT(not modified());
            MADNESS_ASSERT(not doleaves);
            MADNESS_ASSERT(coeff.dim(0)==2*k);
            MADNESS_ASSERT(2*NDIM==coeff.ndim());

            const SeparatedConvolutionData<Q,NDIM>* op = getop(source.level(), shift, source);

            // prepare access to the singular vectors
            std::vector<Slice> s(coeff.config().dim_per_vector()+1,_);
            // can't use predefined slices and vectors -- they have the wrong dimension
            const std::vector<Slice> s00(coeff.ndim(),Slice(0,k-1));

            // some workspace
            Tensor<resultT> work1(v2k,false), work2(v2k,false);
            Tensor<Q> work5(2*k,2*k);

            // sliced input and final result
            const GenTensor<T> f0 = copy(coeff(s00));
            GenTensor<resultT> final=copy(coeff);
            GenTensor<resultT> final0=copy(f0);

            tol = tol/rank*0.01; // Error is per separated term
            tol2= tol2/rank;

            for (int r=0; r<coeff.rank(); ++r) {

                // get the appropriate singular vector (left or right depends on particle)
                // and apply the full tensor muopxv_fast on it, term by term
                s[0]=Slice(r,r);
                const Tensor<T> chunk=coeff.config().ref_vector(particle()-1)(s).reshape(2*k,2*k,2*k);
                const Tensor<T> chunk0=f0.config().ref_vector(particle()-1)(s).reshape(k,k,k);
//                const double weight=coeff.config().weights(r);

                // accumulate all terms of the operator for a specific term of the function
                Tensor<resultT> result(v2k), result0(vk);

                ApplyTerms at;
                at.r_term=true;
                at.t_term=source.level()>0;

                // this loop will return on result and result0 the terms [(P+Q) G (P+Q)]_1,
                // and [P Q P]_1, respectively
                for (int mu=0; mu<rank; ++mu) {
                    const SeparatedConvolutionInternal<Q,NDIM>& muop =  op->muops[mu];

//                    if (muop.norm > tol2*std::abs(weight)) {

                        Q fac = ops[mu].getfac();
                        muopxv_fast(at, muop.ops, chunk, chunk0, result, result0,
                                tol/std::abs(fac), fac, work1, work2, work5);

//                    }
                }


                // reinsert the transformed terms into result, leaving the other particle unchanged
                MADNESS_ASSERT(final.config().has_structure());
                final.config().ref_vector(particle()-1)(s)=result;

                if (source.level()>0) {
                    final0.config().ref_vector(particle()-1)(s)=result0;
                } else {
                    final0.config().ref_vector(0)(s)=0.0;
                    final0.config().ref_vector(1)(s)=0.0;
                }

            }

            final(s00)+=final0;
            final.reduce_rank(tol2);

            return final;
        }



        template <typename T>
        GenTensor<TENSOR_RESULT_TYPE(T,Q)> apply2(const Key<NDIM>& source,
                                              const Key<NDIM>& shift,
                                              const GenTensor<T>& coeff,
                                              double tol, double tol2) const {
            PROFILE_MEMBER_FUNC(SeparatedConvolution);
            typedef TENSOR_RESULT_TYPE(T,Q) resultT;

            MADNESS_ASSERT(coeff.ndim()==NDIM);
            MADNESS_ASSERT(coeff.tensor_type()==TT_2D); // we use the rank below
//            MADNESS_EXCEPTION("no apply2",1);
            const TensorType tt=coeff.tensor_type();

            const GenTensor<T>* input = &coeff;
            GenTensor<T> dummy;

            if (not modified()) {
                if (coeff.dim(0) == k) {
                    // This processes leaf nodes with only scaling
                    // coefficients ... FuncImpl::apply by default does not
                    // apply the operator to these since for smoothing operators
                    // it is not necessary.  It is necessary for operators such
                    // as differentiation and time evolution and will also occur
                    // if the application of the operator widens the tree.
                    dummy = GenTensor<T>(v2k,coeff.tensor_type());
                    dummy(s0) += coeff;
                    input = &dummy;
                }
                else {
                    MADNESS_ASSERT(coeff.dim(0)==2*k);
                }
            }

            tol = tol/rank; // Error is per separated term
            tol2= tol2/rank;

            const SeparatedConvolutionData<Q,NDIM>* op = getop(source.level(), shift, source);

            GenTensor<resultT> r, r0, result, result0;
            GenTensor<resultT> work1(v2k,tt), work2(v2k,tt);
            GenTensor<Q> work5(v2k,tt);

            if (modified()) {
                r=GenTensor<resultT>(vk,tt);
                work1=GenTensor<resultT>(vk,tt);
                work2=GenTensor<resultT>(vk,tt);
                work5=GenTensor<Q>();

            }

            // collect the results of the individual operator terms
            std::list<GenTensor<T> > r_list;
            std::list<GenTensor<T> > r0_list;

//            const GenTensor<T> f0 = copy(coeff(s0));
            const GenTensor<T> f0 = copy((*input)(s0));
            for (int mu=0; mu<rank; ++mu) {
                const SeparatedConvolutionInternal<Q,NDIM>& muop =  op->muops[mu];
                //print("muop",source, shift, mu, muop.norm);

                // delta(g)  <  delta(T) * || f ||
                if (muop.norm > tol) {

                    // get maximum rank of coeff to contribute:
                    //  delta(g)  <  eps  <  || T || * delta(f)
                    //  delta(coeff) * || T || < tol2
                        const int r_max=SRConf<T>::max_sigma(tol2/muop.norm,coeff.rank(),coeff.config().weights_);
                    //                  print("r_max",coeff.config().weights(r_max));

                        // note that max_sigma is inclusive!
                    if (r_max>=0) {
                        const GenTensor<resultT> chunk=input->get_configs(0,r_max);
                        const GenTensor<resultT> chunk0=f0.get_configs(0,r_max);

                        double cpu0=cpu_time();

                        Q fac = ops[mu].getfac();
                        muopxv_fast2(source.level(), muop.ops, chunk, chunk0, r, r0,
                                tol/std::abs(fac), fac, work1, work2, work5);
                        double cpu1=cpu_time();
                        timer_low_transf.accumulate(cpu1-cpu0);

                        r_list.push_back(r);
                        r0_list.push_back(r0);
                    }
                }
            }

            // finally accumulate all the resultant terms into one tensor
            double cpu0=cpu_time();

            result0=reduce(r0_list,tol2*rank);
            if (r_list.size()>0) r_list.front()(s0)+=result0;
            result=reduce(r_list,tol2*rank);
            result.reduce_rank(tol2*rank);

            double cpu1=cpu_time();
            timer_low_accumulate.accumulate(cpu1-cpu0);
            return result;
        }


        template<typename T>
        double estimate_costs(const Key<NDIM>& source,
                const Key<NDIM>& shift,
                const GenTensor<T>& coeff,
                double tol, double tol2) const {

            if (coeff.tensor_type()==TT_FULL) return 0.5;
            if (2*NDIM==coeff.ndim()) return 1.5;
            MADNESS_ASSERT(NDIM==coeff.ndim());
            MADNESS_ASSERT(coeff.tensor_type()==TT_2D);

            const SeparatedConvolutionData<Q,NDIM>* op = getop(source.level(), shift, source);

            tol = tol/rank; // Error is per separated term
            tol2= tol2/rank;

            const double full_operator_cost=pow(coeff.dim(0),NDIM+1);
            const double low_operator_cost=pow(coeff.dim(0),NDIM/2+1);
            const double low_reduction_cost=pow(coeff.dim(0),NDIM/2);

            double full_cost=0.0;
            double low_cost=0.0;

            for (int mu=0; mu<rank; ++mu) {
                const SeparatedConvolutionInternal<Q,NDIM>& muop =  op->muops[mu];

                // delta(g)  <  delta(T) * || f ||
                if (muop.norm > tol) {
                        // note that max_sigma is inclusive: it returns a slice w(Slice(0,i))
                    long nterms=SRConf<T>::max_sigma(tol2/muop.norm,coeff.rank(),coeff.config().weights_)+1;

                    // take only the first overlap computation of rank reduction into account
                    low_cost+=nterms*low_operator_cost + 2.0*nterms*nterms*low_reduction_cost;

                    full_cost+=full_operator_cost;
                }
            }

            // include random empirical factor of 2
            double ratio=-1.0;
            if (low_cost>0.0) ratio=full_cost/low_cost;
//            print("nterms, full, low, full/low", full_cost, low_cost,shift.distsq(), ratio);
            return ratio;

        }

    };



    static
    inline
    SeparatedConvolution<double_complex,3> PeriodicHFExchangeOperator(World& world,
                                                   Vector<double,3> args,
                                                   double lo,
                                                   double eps,
                                                   const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
                                                   int k=FunctionDefaults<3>::get_k())
    {
        const Tensor<double>& cell_width = FunctionDefaults<3>::get_cell_width();
        double hi = cell_width.normf(); // Diagonal width of cell
        if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation

        GFit<double,3> fit=GFit<double,3>::CoulombFit(lo,hi,eps,false);
                Tensor<double> coeff=fit.coeffs();
                Tensor<double> expnt=fit.exponents();

                if (bc(0,0) == BC_PERIODIC) {
                fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), true);
        }

        return SeparatedConvolution<double_complex,3>(world, args, coeff, expnt, bc, k, false);
//        return SeparatedConvolution<double_complex,3>(world, coeff, expnt, bc, k);

    }

    static
    inline
    SeparatedConvolution<double,3> CoulombOperator(World& world,
                                                   double lo,
                                                   double eps,
                                                   const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
                                                   int k=FunctionDefaults<3>::get_k())
    {
        const Tensor<double>& cell_width = FunctionDefaults<3>::get_cell_width();
        double hi = cell_width.normf(); // Diagonal width of cell
        if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation

        GFit<double,3> fit=GFit<double,3>::CoulombFit(lo,hi,eps,false);
                Tensor<double> coeff=fit.coeffs();
                Tensor<double> expnt=fit.exponents();


        if (bc(0,0) == BC_PERIODIC) {
            fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), true);
        }
        return SeparatedConvolution<double,3>(world, coeff, expnt, bc, k);
    }


    static
    inline
    SeparatedConvolution<double,3>* CoulombOperatorPtr(World& world,
                                                       double lo,
                                                       double eps,
                                                       const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
                                                       int k=FunctionDefaults<3>::get_k())
    {
        const Tensor<double>& cell_width = FunctionDefaults<3>::get_cell_width();
        double hi = cell_width.normf(); // Diagonal width of cell
        if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation
        GFit<double,3> fit=GFit<double,3>::CoulombFit(lo,hi,eps,false);
                Tensor<double> coeff=fit.coeffs();
                Tensor<double> expnt=fit.exponents();

        if (bc(0,0) == BC_PERIODIC) {
            fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), true);
        }
        return new SeparatedConvolution<double,3>(world, coeff, expnt, bc, k);
    }


    template <std::size_t NDIM>
    static
    inline
    SeparatedConvolution<double,NDIM> BSHOperator(World& world,
                                                  double mu,
                                                  double lo,
                                                  double eps,
                                                  const BoundaryConditions<NDIM>& bc=FunctionDefaults<NDIM>::get_bc(),
                                                  int k=FunctionDefaults<NDIM>::get_k())
    {
        if (eps>1.e-4) {
                if (world.rank()==0) print("the accuracy in BSHOperator is too small, tighten the threshold",eps);
                MADNESS_EXCEPTION("0",1);
        }
        const Tensor<double>& cell_width = FunctionDefaults<NDIM>::get_cell_width();
        double hi = cell_width.normf(); // Diagonal width of cell
        if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation

        GFit<double,NDIM> fit=GFit<double,NDIM>::BSHFit(mu,lo,hi,eps,false);
                Tensor<double> coeff=fit.coeffs();
                Tensor<double> expnt=fit.exponents();

        if (bc(0,0) == BC_PERIODIC) {
            fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false);
        }

        return SeparatedConvolution<double,NDIM>(world, coeff, expnt, bc, k);
    }


    static
    inline
    SeparatedConvolution<double,3> BSHOperator3D(World& world,
                                                 double mu,
                                                 double lo,
                                                 double eps,
                                                 const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
                                                 int k=FunctionDefaults<3>::get_k())

    {
        const Tensor<double>& cell_width = FunctionDefaults<3>::get_cell_width();
        double hi = cell_width.normf(); // Diagonal width of cell
        if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation

        GFit<double,3> fit=GFit<double,3>::BSHFit(mu,lo,hi,eps,false);
                Tensor<double> coeff=fit.coeffs();
                Tensor<double> expnt=fit.exponents();

        if (bc(0,0) == BC_PERIODIC) {
            fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false);
        }
        return SeparatedConvolution<double,3>(world, coeff, expnt, bc, k);
    }

    static inline SeparatedConvolution<double,3> SlaterF12Operator(World& world,
                double mu, double lo, double eps,
                const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
                int k=FunctionDefaults<3>::get_k()) {

        const Tensor<double>& cell_width = FunctionDefaults<3>::get_cell_width();
        double hi = cell_width.normf(); // Diagonal width of cell
        if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation

        GFit<double,3> fit=GFit<double,3>::SlaterFit(mu,lo,hi,eps,false);
                Tensor<double> coeff=fit.coeffs();
                Tensor<double> expnt=fit.exponents();

        if (bc(0,0) == BC_PERIODIC) {
            fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false);
        }
        return SeparatedConvolution<double,3>(world, coeff, expnt, bc, k, false, mu);
    }


    static
    inline
    SeparatedConvolution<double,3>* BSHOperatorPtr3D(World& world,
                                                     double mu,
                                                     double lo,
                                                     double eps,
                                                     const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
                                                     int k=FunctionDefaults<3>::get_k())
    {
        const Tensor<double>& cell_width = FunctionDefaults<3>::get_cell_width();
        double hi = cell_width.normf(); // Diagonal width of cell
        if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation

        GFit<double,3> fit=GFit<double,3>::BSHFit(mu,lo,hi,eps,false);
                Tensor<double> coeff=fit.coeffs();
                Tensor<double> expnt=fit.exponents();

        if (bc(0,0) == BC_PERIODIC) {
            fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false);
        }
        return new SeparatedConvolution<double,3>(world, coeff, expnt, bc, k);
    }



    static
    inline
    std::vector< std::shared_ptr< SeparatedConvolution<double,3> > >
    GradCoulombOperator(World& world,
                        double lo,
                        double eps,
                        const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
                        int k=FunctionDefaults<3>::get_k())
    {
        typedef SeparatedConvolution<double,3> real_convolution_3d;
        typedef std::shared_ptr<real_convolution_3d> real_convolution_3d_ptr;
        const double pi = constants::pi;
        const Tensor<double> width = FunctionDefaults<3>::get_cell_width();
        double hi = width.normf(); // Diagonal width of cell
        const bool isperiodicsum = (bc(0,0)==BC_PERIODIC);
        if (isperiodicsum) hi *= 100; // Extend range for periodic summation

        GFit<double,3> fit=GFit<double,3>::CoulombFit(lo,hi,eps,false);
                Tensor<double> coeff=fit.coeffs();
                Tensor<double> expnt=fit.exponents();

        if (bc(0,0) == BC_PERIODIC) {
            fit.truncate_periodic_expansion(coeff, expnt, width.max(), true);
        }

        int rank = coeff.dim(0);

        std::vector<real_convolution_3d_ptr> gradG(3);

        for (int dir=0; dir<3; dir++) {
            std::vector< ConvolutionND<double,3> > ops(rank);
            for (int mu=0; mu<rank; mu++) {
                // We cache the normalized operator so the factor is the value we must multiply
                // by to recover the coeff we want.
                double c = std::pow(sqrt(expnt(mu)/pi),3); // Normalization coeff
                ops[mu].setfac(coeff(mu)/c/width[dir]);

                for (int d=0; d<3; d++) {
                    if (d != dir)
                        ops[mu].setop(d,GaussianConvolution1DCache<double>::get(k, expnt(mu)*width[d]*width[d], 0, isperiodicsum));
                }
                ops[mu].setop(dir,GaussianConvolution1DCache<double>::get(k, expnt(mu)*width[dir]*width[dir], 1, isperiodicsum));
            }
            gradG[dir] = real_convolution_3d_ptr(new SeparatedConvolution<double,3>(world, ops));
        }

        return gradG;
    }

    namespace archive {
        template <class Archive, class T, std::size_t NDIM>
        struct ArchiveLoadImpl<Archive,const SeparatedConvolution<T,NDIM>*> {
            static inline void load(const Archive& ar, const SeparatedConvolution<T,NDIM>*& ptr) {
                WorldObject< SeparatedConvolution<T,NDIM> >* p = NULL;
                ar & p;
                ptr = static_cast< const SeparatedConvolution<T,NDIM>* >(p);
            }
        };

        template <class Archive, class T, std::size_t NDIM>
        struct ArchiveStoreImpl<Archive,const SeparatedConvolution<T,NDIM>*> {
            static inline void store(const Archive& ar, const SeparatedConvolution<T,NDIM>*const& ptr) {
                ar & static_cast< const WorldObject< SeparatedConvolution<T,NDIM> >* >(ptr);
            }
        };
    }

}




#endif // MADNESS_MRA_OPERATOR_H__INCLUDED