Interpolate1D "in place" allowed ?

[Guillaume_A_2]

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

[mecej4]

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.

[VictoriyaS_F_Intel · Accepted solution]

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

[Guillaume_A_2]

Ok, Thanks !

[VictoriyaS_F_Intel]

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.