Interpolate1D "in place" allowed ?

Guillaume_A_2 — Mon, 16 Jan 2017 10:50:19 GMT

Hi everyone,

I notice that df?Interpolate1D and df?InterpolateEx1D have different behaviour when working InPlace (&site == &r) : InterpolateEx1D works, Interpolate1D does not.

Here is my sample test:

int flatExtrapolation(long long* n, long long* /*cell*/, double* /*site*/, double* r, void* params)
{
  const double value = *static_cast<double*>(params);
  std::fill(r, r + *n, value);
  return 0;
}
class ExtrapolationFunction
{
public:
  virtual ~ExtrapolationFunction(){}
  explicit ExtrapolationFunction(double value)
  {
    _extrapolationFunction = static_cast<dfdInterpCallBack>(flatExtrapolation);
    _value = value;
  }
  const dfdInterpCallBack& getFunction() const { return _extrapolationFunction; }
  const void* getParameters() const { return &_value; }
private:
  dfdInterpCallBack _extrapolationFunction;
  double _value;
};

void interpolInPlaceTest()
{
  const int lengthX = 7;
  const double minAbscissaRef = 0.0;
  const double maxAbscissaRef = 3.0;
  const int lengthOrdinateRef = 4;

  double *X = new double[lengthX];
  double *XBackUp = new double[lengthX];
  double *Z = new double[lengthX];

  X[0] = 0.0;
  X[1] = 0.5;
  X[2] = 1.0;
  X[3] = 1.5;
  X[4] = 2.0;
  X[5] = 2.5;
  X[6] = 3.0;
  for (int i = 0; i < lengthX; i++)
  {
    XBackUp = X;
    Z = 0.0;
  }
  
  double* ordinates = new double[lengthOrdinateRef];
  double* abscissae = new double[lengthOrdinateRef];
  double step = (maxAbscissaRef - minAbscissaRef) / (lengthOrdinateRef - 1);
  for (int i = 0; i < lengthOrdinateRef; i++)
  {
    ordinates = i + 20.0;
    abscissae = minAbscissaRef + step*i;
  }

  MKL_INT ny = 1;
  DFTaskPtr task = nullptr;
  double* splineCoeffs = new double[DF_PP_LINEAR*(lengthOrdinateRef - 1)];

  dfdNewTask1D(&task, lengthOrdinateRef, abscissae, DF_NON_UNIFORM_PARTITION, ny, ordinates, DF_MATRIX_STORAGE_ROWS);
  dfdEditPPSpline1D(task, DF_PP_LINEAR, DF_PP_DEFAULT, DF_NO_BC, nullptr, DF_NO_IC, nullptr, splineCoeffs, DF_NO_HINT);
  dfdConstruct1D(task, DF_PP_SPLINE, DF_METHOD_STD);

  int ndorder = 1;
  MKL_INT* dorder = new MKL_INT[ndorder];
  dorder[0] = 1;

  // ****** Interpolation without extrapolation ******
  int status = dfdInterpolate1D(task, DF_INTERP, DF_METHOD_PP, lengthX, X, DF_NON_UNIFORM_PARTITION, ndorder, dorder,
    nullptr, Z, DF_MATRIX_STORAGE_ROWS, nullptr);
  printf("\ndfdInterpolate1D Out of Place status: %d\n", status);
  for (int i = 0; i < lengthX; i++)
    printf("%d\t%lf\t%lf\n", i, XBackUp, Z);

  status = dfdInterpolate1D(task, DF_INTERP, DF_METHOD_PP, lengthX, X, DF_NON_UNIFORM_PARTITION, ndorder, dorder,
    nullptr, X, DF_MATRIX_STORAGE_ROWS, nullptr);
  printf("\ndfdInterpolate1D In Place status: %d\n", status);
  for (int i = 0; i < lengthX; i++)
    printf("%d\t%lf\t%lf\n", i, XBackUp, X);

  // ****** Interpolation with extrapolation ******
  // reset buffers
  for (int i = 0; i < lengthX; i++)
  {
    X = XBackUp;
    Z = 0.0;
  }
  ExtrapolationFunction* _extrapolationFunc = new ExtrapolationFunction(42.0);
  status = dfdInterpolateEx1D(task, DF_INTERP, DF_METHOD_PP, lengthX, X, DF_NON_UNIFORM_PARTITION, ndorder, dorder,
    nullptr, Z, DF_MATRIX_STORAGE_ROWS, nullptr,
    _extrapolationFunc->getFunction(),
    _extrapolationFunc->getParameters(), 
    _extrapolationFunc->getFunction(),
    _extrapolationFunc->getParameters(), 
    nullptr, nullptr, // interpolation function 
    nullptr, nullptr);
  printf("\ndfdInterpolateEx1D Out of Place status: %d\n", status);
  for (int i = 0; i < lengthX; i++)
    printf("%d\t%lf\t%lf\n", i, XBackUp, Z);

  status = dfdInterpolateEx1D(task, DF_INTERP, DF_METHOD_PP, lengthX, X, DF_NON_UNIFORM_PARTITION, ndorder, dorder,
    nullptr, X, DF_MATRIX_STORAGE_ROWS, nullptr,
    _extrapolationFunc->getFunction(),
    _extrapolationFunc->getParameters(), 
    _extrapolationFunc->getFunction(),
    _extrapolationFunc->getParameters(),
    nullptr, nullptr, // interpolation function 
    nullptr, nullptr);
  printf("\ndfdInterpolateEx1D In Place status: %d\n", status);
  for (int i = 0; i < lengthX; i++)
    printf("%d\t%lf\t%lf\n", i, XBackUp, X);

  delete[] X;
  delete[] XBackUp;
  delete[] Z;
  delete[] ordinates;
  delete[] abscissae;
  delete[] splineCoeffs;
  delete[] dorder;
  delete _extrapolationFunc;
  dfDeleteTask(&task);
}

And here is the output :

dfdInterpolate1D Out of Place status: 0
0   0.000000   20.000000
1   0.500000   20.500000
2   1.000000   21.000000
3   1.500000   21.500000
4   2.000000   22.000000
5   2.500000   22.500000
6   3.000000   23.000000
dfdInterpolate1D In Place status: 0
0   0.000000   20.000000
1   0.500000   20.000000
2   1.000000   20.000000
3   1.500000   20.000000
4   2.000000   20.000000
5   2.500000   20.000000
6   3.000000   20.000000
dfdInterpolateEx1D Out of Place status: 0
0   0.000000   20.000000
1   0.500000   20.500000
2   1.000000   21.000000
3   1.500000   21.500000
4   2.000000   22.000000
5   2.500000   22.500000
6   3.000000   42.000000
dfdInterpolateEx1D In Place status: 0
0   0.000000   20.000000
1   0.500000   20.500000
2   1.000000   21.000000
3   1.500000   21.500000
4   2.000000   22.000000
5   2.500000   22.500000
6   3.000000   42.000000

When using dfdInterpolate1D In Place the status is OK but the results are wrong. That were the results with Composer 2016 (first release).

If I run the same code with Composer 2017 Update 1, the output of dfdInterpolateEx1D In Place become wrong too !

So, I am a bit lost... Is working In Place with dfdInterpolate1D (or dfdInterpolateEx1D) allowed or not ? Is that supposed to work ? I did not find the answer in the documentation.

Thanks,

Guix

In general, it is not correct

mecej4 — Mon, 16 Jan 2017 12:49:43 GMT

In general, it is not correct to alias one subroutine argument to another argument, unless the documentation explicitly lists arguments that may be aliased.

The Fortran standard prohibits such aliasing of arguments, and being able to rely on arguments being distinct is what enables many code optimizations to be performed by the compiler. If you are calling a library routine that was compiled from Fortran sources, the rule applies to calls made to the routine from code written in any language.

Hello Guix,

VictoriyaS_F_Intel — Tue, 17 Jan 2017 08:20:25 GMT

Hello Guix,

The functions of Data Fitting component do not support in-place computations.

It is a luck that you get the correct results with in-place invocation of dfdInterpolateEx1D. Please use non-aliasing memory locations for input and output parameters of Data Fitting functions.

Best regards,

Victoriya

Ok, Thanks !

Guillaume_A_2 — Tue, 17 Jan 2017 08:48:39 GMT

Ok, Thanks !

We plan to add a note into

VictoriyaS_F_Intel — Tue, 17 Jan 2017 12:47:22 GMT

We plan to add a note into Data Fitting documentation. The note will explicitly specify that the in-place usage of df?Interpolate1D/df?InterpolateEx1D is not supported.

topic In general, it is not correct in Intel® oneAPI Math Kernel Library

Interpolate1D "in place" allowed ?

In general, it is not correct

Hello Guix,

Ok, Thanks !

We plan to add a note into