Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[oneMKL] add Data Fitting domain to oneMKL #413

Closed
wants to merge 10 commits into from
Closed
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 12 additions & 11 deletions source/elements/oneMKL/source/architecture/api_design.inc.rst
Original file line number Diff line number Diff line change
Expand Up @@ -16,17 +16,18 @@ oneMKL namespaces

The oneMKL library uses C++ namespaces to organize routines by mathematical domain. All oneMKL objects and routines shall be contained within the ``oneapi::mkl`` base namespace. The individual oneMKL domains use a secondary namespace layer as follows:

======================== =======================================================================================================
namespace oneMKL domain or content
======================== =======================================================================================================
``oneapi::mkl`` oneMKL base namespace, contains general oneMKL data types, objects, exceptions and routines
``oneapi::mkl::blas`` Dense linear algebra routines from BLAS and BLAS like extensions. The oneapi::mkl::blas namespace should contain two namespaces column_major and row_major to support both matrix layouts. See :ref:`onemkl_blas`
``oneapi::mkl::lapack`` Dense linear algebra routines from LAPACK and LAPACK like extensions. See :ref:`onemkl_lapack`
``oneapi::mkl::sparse`` Sparse linear algebra routines from Sparse BLAS and Sparse Solvers. See :ref:`onemkl_sparse_linear_algebra`
``oneapi::mkl::dft`` Discrete and fast Fourier transformations. See :ref:`onemkl_dft`
``oneapi::mkl::rng`` Random number generator routines. See :ref:`onemkl_rng`
``oneapi::mkl::vm`` Vector mathematics routines, e.g. trigonometric, exponential functions acting on elements of a vector. See :ref:`onemkl_vm`
======================== =======================================================================================================
=========================================== =========================================================================
namespace oneMKL domain or content
=========================================== =========================================================================
``oneapi::mkl`` oneMKL base namespace, contains general oneMKL data types, objects, exceptions and routines
``oneapi::mkl::blas`` Dense linear algebra routines from BLAS and BLAS like extensions. The oneapi::mkl::blas namespace should contain two namespaces column_major and row_major to support both matrix layouts. See :ref:`onemkl_blas`
``oneapi::mkl::lapack`` Dense linear algebra routines from LAPACK and LAPACK like extensions. See :ref:`onemkl_lapack`
``oneapi::mkl::sparse`` Sparse linear algebra routines from Sparse BLAS and Sparse Solvers. See :ref:`onemkl_sparse_linear_algebra`
``oneapi::mkl::dft`` Discrete and fast Fourier transformations. See :ref:`onemkl_dft`
``oneapi::mkl::rng`` Random number generator routines. See :ref:`onemkl_rng`
``oneapi::mkl::vm`` Vector mathematics routines, e.g. trigonometric, exponential functions acting on elements of a vector. See :ref:`onemkl_vm`
``oneapi::mkl::experimental::data_fitting`` Data fitting routines, e.g. interpolate. See :ref:`data_fitting`
andreyfe1 marked this conversation as resolved.
Show resolved Hide resolved
=========================================== =========================================================================

.. note::
:name: Implementation Requirement
Expand Down
90 changes: 90 additions & 0 deletions source/elements/oneMKL/source/domains/data_fitting/cubic.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,90 @@
.. _cubic:

Cubic Splines
=============

.. contents::
:local:
:depth: 2

Cubic splines are splines whose degree is equal to 3.

Cubic splines are described by the following polynomial

.. math::
P_i\left( x \right) = c_{1,i}+ c_{2,i}\left( x - x_i \right) + c_{3,i}{\left( x - x_i \right)}^2+ c_{4,i}{\left( x - x_i \right)}^3,

where

.. math::
x \in \left[ x_i, x_{i+1} \right),

.. math::
i = 1,\cdots , n-1.

There are a lot of different types of cubic splines: Hermite, natural, Akima, Bessel.
andreyfe1 marked this conversation as resolved.
Show resolved Hide resolved
However, the current version of DPC++ API supports only one type: Hermite.

Header File
-----------

.. code:: cpp

#include<oneapi/mkl/experimental/data_fitting.hpp>

Namespace
---------

.. code:: cpp

oneapi::mkl::experimental::data_fitiing

Hermite Spline
--------------

Coefficients of Hermite spline are calculated using the following formulas:
andreyfe1 marked this conversation as resolved.
Show resolved Hide resolved

.. math::
c_{1,i} = f\left( x_i \right),

.. math::
c_{2,i} = s_i,

.. math::
c_{3,i} = \left( \left[ x_i, x_{i+1} \right]f - s_i \right) / \left( \Delta x_i \right) - c_{4,i}\left( \Delta x_i \right),

.. math::
c_{4,i} = \left( s_i + s_{i+1} - 2\left[ x_i, x_{i+1} \right]f \right) / {\left( \Delta x_i \right)}^2,

.. math::
s_i = f^{\left( 1 \right)}\left( x_i \right).

The following boundary conditions are supported for Hermite spline:

- Free end (:math:`f^{(2)}(x_1) = f^{(2)}(x_n) = 0`).
- Periodic.
- First derivative.
- Second Derivative.

Syntax
^^^^^^

.. code:: cpp

namespace cubic_spline {
struct hermite {};
}

Example
^^^^^^^

To create a cubic Hermite spline object use the following:

.. code:: cpp

spline<float, cubic_spline::hermite> val(
/*SYCL queue object*/q,
/*number of spline functions*/ny
);

Follow the :ref:`examples` section to see more complicated examples.
Original file line number Diff line number Diff line change
@@ -0,0 +1,30 @@
.. _data-fitting:

Data Fitting
============

|IONE-MKL| provides spline-based interpolation capabilities that can be used for
spline construction (Linear, Cubic, Quadratic etc.),
to perform cell-search operations, and to approximate functions,
function derivatives, or integrals.

APIs are experimental. It means that no API or ABI backward compatibility are guaranteed.

APIs are based on SYCL USM (Unfied Shared Memory) input data.

Routines
--------

Splines:
:ref:`splines`
Interpolate function:
:ref:`interpolate`

.. toctree::
:maxdepth: 1
:hidden:

data_fitting/terms
data_fitting/splines
data_fitting/interpolate
data_fitting/examples
10 changes: 10 additions & 0 deletions source/elements/oneMKL/source/domains/data_fitting/examples.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
.. _examples:

Examples
========

The following example demonstrates how to construct the linear spline and perform the interpolation.

.. literalinclude:: /_examples/data_fitting_linear.cpp
:language: cpp
:linenos:
106 changes: 106 additions & 0 deletions source/elements/oneMKL/source/domains/data_fitting/interpolate.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,106 @@
.. _interpolate:

Interpolate Function
====================

.. contents::
:local:
:depth: 1

Interpolate function performs computations of function and derivatives values at interpolation sites.

If the sites do not belong to interpolation interval ``[a, b]``, the library uses:

- interpolant :math:`I_0` coefficients computed for interval :math:`[x_0, x_1)` for the
computations at the sites to the left of ``a``.
- interpolant :math:`I_{n-2}` coefficients computed for interval
:math:`[x_{n-2}, x_{n-1})` for the computations at the sites to the right of ``b``.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You mean to say that you do allow for extrapolation, but it uses the interpolant function defined on the closest sub-interval in the interpolation interval :math:[a,b)?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This seems different than your definition of site in the terms.rst

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You mean to say that you do allow for extrapolation, but it uses the interpolant function defined on the closest sub-interval in the interpolation interval :math:[a,b)?

Correct.

This seems different than your definition of site in the terms.rst

Statements here are correct. Will rework "terms.rst"


Interpolation algorithm depends on interpolant's type (e.g., for cubic spline
andreyfe1 marked this conversation as resolved.
Show resolved Hide resolved
interpoilation evaluation of third-order polynomial is performed to obtain function values).
andreyfe1 marked this conversation as resolved.
Show resolved Hide resolved

Header File
-----------

.. code:: cpp

#include<oneapi/mkl/experimental/data_fitting.hpp>

Namespace
---------

.. code:: cpp

oneapi::mkl::experimental::data_fitiing
andreyfe1 marked this conversation as resolved.
Show resolved Hide resolved

Syntax
------

.. code:: cpp

template <typename Interpolant>
sycl::event interpolate(
Interpolant& interpolant,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

are these missing the queue as argument?

Also, do you fit into non-member or member function cases as seen in https://github.com/oneapi-src/oneAPI-spec/blob/main/source/elements/oneMKL/source/architecture/execution_model.inc.rst#member-functions ? If you are of member-function approach, you may want to add some language there to include Data fitting in list with DFt and RNG.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

queue is not missed it's inside interpolant (already associated with interpolant). We also have interfaces with "queue" below.

It's a non-member function. Spec says:

Each oneMKL non-member computational routine takes a sycl::queue reference as its first parameter

I see the point. We need to say somehow that the first argument contains execution context.
Given our case maybe we can rephrase the sentence from execution_model.inc.rst:

Each oneMKL non-member computational routine takes a value that contain execution context as its first parameter

What do you think?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What about changing the spec to say something like this:

In general, each oneMKL non-member computational routine takes a sycl::queue reference as its first parameter.
In some exceptional cases, like in Data Fitting, it may instead take an object as first parameter that itself contains a reference to the sycl execution context.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why is a queue not passed in? why is a context sufficient ? Is there a future call which adds the queue ? How does work get done ?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually the queue is inside (See here: https://github.com/oneapi-src/oneAPI-spec/pull/413/files/#diff-2361f70703f03fec4684fd287b54dae6e0b303c2f095bab62be9b5bc536b51f8R34-R41).
Sorry, I might misspell it trying to explain it more generally.

typename Interpolant::fp_type* sites,
std::int64_t n_sites,
typename Interpolant::fp_type* results,
const std::vector<sycl::event>& dependencies,
interpolate_hint ResultHint = interpolate_hint::funcs_sites_ders,
site_hint SiteHint = site_hint::non_uniform); // (1)

template <typename Interpolant>
sycl::event interpolate(
Interpolant& interpolant,
typename Interpolant::fp_type* sites,
std::int64_t n_sites,
typename Interpolant::fp_type* results,
std::bitset<32> der_indicator,
Copy link
Contributor

@spencerpatty spencerpatty Mar 17, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is there any reason you shortened this from derivative_indicator to der_indicator ? The former would make more sense to me since almost everything else is already full named.

Copy link
Contributor

@spencerpatty spencerpatty Mar 17, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

although I suppose you have used "ders" in the interpolate_hint types... so maybe there is some precedence.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You're right that we use "ders" as a part of interpolation hint names. Am I right that you're ok if we use ders_indicator as a name?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I suppose so, but if I were you, I would choose a general pattern so it is predictable for all names... lower snake case with full words ... I see you sometimes shortening words sometimes not shortening them... is there a pattern to it? maybe you wan to rethink all the names into a single pattern :)

See if you can write down a general naming pattern: recall this spec is supposed to be written so that another developer could implement the library, so you want to describe naming patterns when applicable in case someone wants to suggest an extension.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

on the other hand, I am ok with just a few shortened cases as long as it is readable... but you should review all your names and see if there is a general pattern that you can (and want to) write down... you may discover you want to be more consistent, or that you are ok with it :) To me this is just one aspect of the full review.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok. I see. Yeah, we need to think about general pattern and exceptions for it (if any). As for me interpolation sites will look too long... However, we definitely need to think about it. Thanks

const std::vector<sycl::event>& dependencies = {},
interpolate_hint ResultHint = interpolate_hint::funcs_sites_ders,
site_hint SiteHint = site_hint::non_uniform); // (2)

template <typename Interpolant>
sycl::event interpolate(
sycl::queue& q,
andreyfe1 marked this conversation as resolved.
Show resolved Hide resolved
const Interpolant& interpolant,
typename Interpolant::fp_type* sites,
std::int64_t n_sites,
typename Interpolant::fp_type* results,
const std::vector<sycl::event>& dependencies,
interpolate_hint ResultHint = interpolate_hint::funcs_sites_ders,
site_hint SiteHint = site_hint::non_uniform); // (3)

template <typename Interpolant>
sycl::event interpolate(
sycl::queue& q,
const Interpolant& interpolant,
typename Interpolant::fp_type* sites,
std::int64_t n_sites,
typename Interpolant::fp_type* results,
std::bitset<32> der_indicator,
const std::vector<sycl::event>& dependencies = {},
interpolate_hint ResultHint = interpolate_hint::funcs_sites_ders,
site_hint SiteHint = site_hint::non_uniform); // (4)

For all functions users can provide ``SiteHint`` and ``ResultHint`` to specify
the layout of ``sites`` and ``results`` respectively.
If ``results`` layout doesn't satisfy ``ResultHint`` and/or
``sites`` layout doesn't satisfy ``SiteHint``, behavior is undefined.
Returns the SYCL event of the submitted task.

#. Performs computations of function values only using the SYCL queue
associated with ``interpolant``.
#. Performs computations of certain derivatives
(function values is considered as a zero derivative) which are indicated in
``der_indicator`` (each bit corresponds to certain derivative starting from lower bit)
using the SYCL queue associated with ``interpolant``.
#. Performs computations of function values only using ``q`` as an input argument
that should be created from the same context and device as the SYCL queue
associated with ``interpolant``.
#. Performs computations of certain derivatives
(function values is considered as a zero derivative) which are indicated in
``der_indicator`` (each bit corresponds to certain derivative starting from lower bit)
using ``q`` as an input argument that should be created from
the same context and device as the SYCL queue associated with ``interpolant``.

Follow the :ref:`examples` section to see examples of the interpolation function usage.
66 changes: 66 additions & 0 deletions source/elements/oneMKL/source/domains/data_fitting/linear.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,66 @@
.. _linear:

Linear Spline
=============

.. contents::
:local:
:depth: 1

Linear spline is a spline whose degree is equal to 1.

It's described by the following polynomial

.. math::
P_i\left( x \right) = c_{1,i} + c_{2,i}\left( x - x_i \right),

where

.. math::
x \in \left[ x_i, x_{i+1} \right),

.. math::
c_{1,i} = f\left( x_i \right),

.. math::
c_{2,i} = \left[ x_i, x_{i+1} \right]f,

.. math::
i = 1, \cdots, n-1.

Header File
-----------

.. code:: cpp

#include<oneapi/mkl/experimental/data_fitting.hpp>

Namespace
---------

.. code:: cpp

oneapi::mkl::experimental::data_fitiing

Syntax
------

.. code:: cpp

namespace linear_spline {
struct default_type {};
}

Example
-------

To create a linear spline object use the following:

.. code:: cpp

spline<float, linear_spline::default_type> val(
/*SYCL queue object*/q,
/*number of spline functions*/ny
);

Follow the :ref:`examples` section to see more complicated examples.
Loading