function template parametrized by other function with different number of arguments

Question

I'm able to make function template parametrized by an other function, however, I don't know how to do it when I want to parametrize it by function with different number of arguments.

See this code:

#include <stdio.h>
#include <math.h>

template < double FUNC( double a ) >
void seq_op( int n, double * as ){
    for (int i=0; i<n; i++){  printf( " %f \n", FUNC( as[i] )  ); }
} 

template < double FUNC( double a, double b ) >
void seq_op_2( int n, double * as, double * bs ){
    for (int i=0; i<n; i++){  printf( " %f \n", FUNC( as[i], bs[i] )  ); }
} 

double a_plus_1  ( double a ){ return a + 1.0; }
double a_sq      ( double a ){ return a*a;     }

double a_plus_b ( double a, double b ){ return a + b; }
double a_times_b( double a, double b ){ return a * b; }


double as[5] = {1,2,3,4};
double bs[5] = {2,2,2,2};

// FUNCTION ======  main
int main(){
    printf( "seq_op   <a_plus_1>  ( 5, as );\n");      seq_op   <a_plus_1>  ( 4, as );
    printf( "seq_op   <a_sq>      ( 5, as );\n");      seq_op   <a_sq>      ( 4, as );
    printf( "seq_op_2 <a_plus_b>  ( 5, as, bs );\n");  seq_op_2 <a_plus_b>  ( 4, as, bs );
    printf( "seq_op_2 <a_times_b> ( 5, as, bs );\n");  seq_op_2 <a_times_b> ( 4, as, bs );
}

is there a way how to make common template for both cases?

Why I need such silly thing? A more practical example are this two functions which differs only in one line:

#define i3D( ix, iy, iz )  ( iz*nxy + iy*nx + ix  ) 

void getLenardJonesFF( int natom, double * Rs_, double * C6, double * C12 ){
    Vec3d * Rs = (Vec3d*) Rs_;
    int nx  = FF::n.x;
    int ny  = FF::n.y;
    int nz  = FF::n.z;
    int nxy = ny * nx;
    Vec3d rProbe;  rProbe.set( 0.0, 0.0, 0.0 ); // we may shift here
    for ( int ia=0; ia<nx; ia++ ){ 
        printf( " ia %i \n", ia );
        rProbe.add( FF::dCell.a );  
        for ( int ib=0; ib<ny; ib++ ){ 
            rProbe.add( FF::dCell.b );
            for ( int ic=0; ic<nz; ic++ ){
                rProbe.add( FF::dCell.c );
                Vec3d f; f.set(0.0,0.0,0.0);
                for(int iatom=0; iatom<natom; iatom++){
                    // only this line differs
                    f.add( forceLJ( Rs[iatom] - rProbe, C6[iatom], C12[iatom] ) );
                }
                FF::grid[ i3D( ia, ib, ic ) ].add( f );
            } 
            rProbe.add_mul( FF::dCell.c, -nz );
        } 
        rProbe.add_mul( FF::dCell.b, -ny );
    }
}

void getCoulombFF( int natom, double * Rs_, double * kQQs ){
    Vec3d * Rs = (Vec3d*) Rs_;
    int nx  = FF::n.x;
    int ny  = FF::n.y;
    int nz  = FF::n.z;
    int nxy = ny * nx;
    Vec3d rProbe;  rProbe.set( 0.0, 0.0, 0.0 ); // we may shift here
    for ( int ia=0; ia<nx; ia++ ){ 
        printf( " ia %i \n", ia );
        rProbe.add( FF::dCell.a );  
        for ( int ib=0; ib<ny; ib++ ){ 
            rProbe.add( FF::dCell.b );
            for ( int ic=0; ic<nz; ic++ ){
                rProbe.add( FF::dCell.c );
                Vec3d f; f.set(0.0,0.0,0.0);
                for(int iatom=0; iatom<natom; iatom++){
                    // only this line differs
                    f.add( forceCoulomb( Rs[iatom] - rProbe, kQQs[iatom] );
                }
                FF::grid[ i3D( ia, ib, ic ) ].add( f );
            } 
            rProbe.add_mul( FF::dCell.c, -nz );
        } 
        rProbe.add_mul( FF::dCell.b, -ny );
    }
}

Answer 1

You should be able to combine the two functions using a combination of std::bind() and std::function() (see code on coliru):

#include <stdio.h>
#include <functional>

using namespace std::placeholders;


double getLJForceAtoms (int, int, double*, double*, double*)
{
    printf("getLJForceAtoms\n");
    return 0;
}


double getCoulombForceAtoms (int, int, double*, double*)
{
    printf("getCoulombForceAtoms\n");
    return 0;
}


void getFF (int natom, double* Rs_, std::function<double(int, int, double*)> GetForce)
{
    int rProbe = 1;

    double Force = GetForce(rProbe, natom, Rs_);
}


int main ()
{
    double* C6 = nullptr;
    double* C12 = nullptr;
    double *kQQs = nullptr;
    double* Rs_ = nullptr;

    auto getLJForceFunc = std::bind(getLJForceAtoms, _1, _2, _3, C6, C12);
    auto getCoulombForceFunc = std::bind(getCoulombForceAtoms, _1, _2, _3, kQQs);

    getFF(1, Rs_, getLJForceFunc);
    getFF(1, Rs_, getCoulombForceFunc);

    return 0;
}

which outputs the expected:

getLJForceAtoms
getCoulombForceAtoms

Update -- On Performance

While it is natural to be concerned about performance of using std::function vs templates I would not omit a possible solution without first benchmarking and profiling it.

I can't compare the performance directly as I would need both your complete source code as well as input data set to make accurate benchmarks but I can do a very simple test to show you what it could look like. If I make the force functions do a little work:

double getLJForceAtoms (int x, int y, double* r1, double* r2, double* r3)
{
    return cos(log2(abs(sin(log(pow(x, 2) + pow(y, 2))))));
}

and then have a very simple getFF() function call them 10 million times I can get a rough comparison between the various design methods (tests done on VS2013, release build, fast optimization flags):

• Direct Call = 1900 ms
• Switch = 1900 ms
• If (flag) = 1900 ms
• Virtual Function = 2400 ms
• std::function = 2400 ms

So the std::function method is about 25% slower in this case but the switch and if methods are the same speed as the direct call case. Depending on how much work your actual force functions do you may get worse or better results. These days, the compiler optimizer and the CPU branch predictor are good enough to do a lot of things that may be surprising or even counter-intuitive, which is why actual testing must be done.

I would do a similar benchmark test with your exact code and data set and see what difference, if any, the various designs have. If you really only have two cases as shown in your question then the "if (flag)" method may be a good choice.