NIntegrate fails when integrating over a list from an external (MathLink) function

Question

NIntegrate fails when integrating over a list from an external (MathLink) function. For simplicity, consider an external function f[x] that returns the list {x,2x}. In Mathematica, the function would be defined as f[x_]={x,2x}. In reality I have a much more complicated function that cannot be implemented inside Mathematica. NIntegrate works fine if the external function f[x] returns a single number, not a list.

Install["f"]
(* LinkObject[...] *)
f[.5]
(* {0.5, 1.} *)
NIntegrate[f[x], {x, 0, 1}]  
(* Integrand f[x] is not numerical at {x} = {0.00795732}. >>*)
(* Integrand f[x] is not numerical at {x} = {0.00795732}. >>*)

This is the file f.tm that defines the external function f[x]:

#include "mathlink.h"

:Begin:
:Function:       f
:Pattern:        f[x_?NumericQ]
:Arguments:      {x}
:ArgumentTypes:  {Real}
:ReturnType:     Manual
:End:

void f(double x)
{
    double y[2];
    y[0] = x;
    y[1] = 2 * x;

    MLPutRealList(stdlink, y, 2);
}

int main(int argc, char* argv[])
{
    return MLMain(argc, argv);
}

Answer 1

NIntegrate can handle vector integrands if they are explicitly in vector form:

f[x_] := {x, 2 x}
NIntegrate[f[x], {x, 0, 1}]
(* {0.5, 1.} *)

When the function is a black box which will only evaluate for numeric arguments, NIntegrate no longer sees the list structure and complains about not getting a scalar value:

g[x_?NumericQ] := {x, 2 x}
NIntegrate[g[x], {x, 0, 1}]
(* NIntegrate::inum: Integrand g[x] is not numerical ... *)

We can still integrate g by supplying the components separately:

g1[x_?NumericQ] := g[x][[1]]
g2[x_?NumericQ] := g[x][[2]]
NIntegrate[{g1[x], g2[x]}, {x, 0, 1}]
(* {0.5, 1.} *)

Answer 2

f[x] should return a single value not a pair. Your example data is of the form

data = {#, 2 #} & /@ Range[0, 1, .05];

Clear[f]

f = Interpolation[data];

Note that f[x] is single-valued

f[.5]

1.

Using NIntegrate

NIntegrate[f[x], {x, 0, 1}]

1.

However, you can Integrate an InterpolatingFunction

Integrate[f[x], {x, 0, 1}]

1.

NIntegrate[f[x], {x, 0, 1}] ==
 Integrate[f[x], {x, 0, 1}] ==
 Integrate[2 x, {x, 0, 1}]

True