DPS921/Intel Math Kernel Library

GPU621/DPS921 | Participants | Groups and Projects | Resources | Glossary

Intel Math Kernel Library

Team Name

Slightly Above Average

Group Members

- Jordan Sie

- James Inkster

- Varshan Nagarajan

Progress

100%/100%

Intel Math Library

Features

Linear Algebra

Intel has implemented linear algebra functions that follow industry standards ([BLAS](http://www.netlib.org/blas/) and [LAPACK](http://www.netlib.org/lapack/)). These functions include those that can do the following:

Level 1: Vector-vector operations

Level 2: Matrix-vector operations

Level 3: Matrix-matrix operations

Sparse Linear Algebra Functions

Able to perform low-level inspector-executor routines on sparse matrices, such as:

- Multiply sparse matrix with dense vector

- Multiply sparse matrix with dense matrix

- Solve linear systems with triangular sparse matrices

- Solve linear systems with general sparse matrices

A sparse matrix is matrix that is mostly empty, these are common in machine learning applications. Using standard linear algebra functions would lead to poor performance and would require greater amounts of storage. Specially written sparse linear algebra functions have better performance and can better compress matrices to save space.

Fast Fourier Transforms

Enabling technology today such as most digital communications, audio compression, image compression, satellite tv, FFT is at the heart of it. A fast Fourier transform (FFT) is an algorithm that computes the discrete Fourier transform (DFT) of a sequence, or its inverse (IDFT).

Random Number Generator

The Intel Math Kernel Library has an interface for RNG routines that use pseudorandom, quasi-random, and non-deterministic generators. These routines are developed the calls to the highly optimized Basic Random Number Generators (BRNGs). All BRNGs differentiate in speeds and properties so its easy to find a optimized one for your application.

All RNG routines can be categorized in several different categories.

- Engines - hold the state of a generator

- Transformation classes - holds different types of statistical distribution

- Generate function - the routine that obtains the random number from the statistical distribution

- Services - using routines that can modify the state of the engine

The generation of numbers is done in 2 steps:

1. generate the state using the engine.

2. iterate over the values and the output is the random numbers.

Data Fitting

The Intel Math Kernal Library provide spline-based interpolation that can be ultilized to approximate functions for derivatives, integrals and cell search operations.

Data Fitting routines use the following workflow to process a task:

- Create a task or multiple tasks.

- Modify the task parameters.

- Perform a Data Fitting computation.

- Destroy the task or tasks.

Data Fitting functions:

- Task Creation and Initialization Routines.

- Task Configuration Routines.

- Computational Routines.

- Task Destructors.

Summary Statistics

The Intel Math Kernal Library has an interface for Summary Statistics that can compute estimates for single, double and multi-dimensional datasets. For example, such parameters may be precision, dimensions of user data, the matrix of the observations, or shapes of data arrays. First you create and intialize the object for the dataset, then you call the summary statistics routine to calculate the estimate.

Summary Statistics calculate:

- Raw and central sums/moments up to the fourth order.

- Variation coefficient.

- Skewness and excess kurtosis.

- Minimum and maximum.

Additional Features:

- Detect outliers in datasets.

- Support missing values in datasets.

- Parameterize correlation matrices.

- Compute quantiles for streaming data.

Vector Math

There are two main set of functions for the Vector Math library that the intel MKL uses they are:

- VM Mathematical Functions�Which allows for it to compute values of mathematical functions e.g. sine, cosine, exponential, or logarithm on vectors that are stored in contiguous memory.

- VM Service Functions are used for showing when catching errors made in the calculations or accuracy. Such as catching error codes or error messages from improper calculations.

Code Samples

These samples are directly from the Intel Math Kernal Library code examples.

Vector Add

``` //============================================================== // Vector Add is the equivalent of a Hello, World! sample for data parallel // programs. Building and running the sample verifies that your development // environment is setup correctly and demonstrates the use of the core features // of DPC++. This sample runs on both CPU and GPU (or FPGA). When run, it // computes on both the CPU and offload device, then compares results. If the // code executes on both CPU and offload device, the device name and a success // message are displayed. And, your development environment is setup correctly! // // For comprehensive instructions regarding DPC++ Programming, go to // https://software.intel.com/en-us/oneapi-programming-guide and search based on // relevant terms noted in the comments. // // DPC++ material used in the code sample: // • A one dimensional array of data shared between CPU and offload device. // • A device queue and kernel. //============================================================== // Copyright © Intel Corporation // // SPDX-License-Identifier: MIT // =============================================================

1. include <CL/sycl.hpp>
2. include <array>
3. include <iostream>
4. if FPGA || FPGA_EMULATOR
5. include <CL/sycl/INTEL/fpga_extensions.hpp>
6. endif

using namespace sycl;

// Array size for this example. constexpr size_t array_size = 10000;

// Create an exception handler for asynchronous SYCL exceptions static auto exception_handler = [](sycl::exception_list e_list) {

 for (std::exception_ptr const &e : e_list) {
   try {
     std::rethrow_exception(e);
   }
   catch (std::exception const &e) {

1. if _DEBUG

     std::cout << "Failure" << std::endl;

1. endif

     std::terminate();
   }
 }

};

//************************************ // Vector add in DPC++ on device: returns sum in 4th parameter "sum". //************************************ void VectorAdd(queue &q, const int *a, const int *b, int *sum, size_t size) {

 // Create the range object for the arrays.
 range<1> num_items{size};

 // Use parallel_for to run vector addition in parallel on device. This
 // executes the kernel.
 //    1st parameter is the number of work items.
 //    2nd parameter is the kernel, a lambda that specifies what to do per
 //    work item. the parameter of the lambda is the work item id.
 // DPC++ supports unnamed lambda kernel by default.
 auto e = q.parallel_for(num_items, [=](auto i) { sum[i] = a[i] + b[i]; });

 // q.parallel_for() is an asynchronous call. DPC++ runtime enqueues and runs
 // the kernel asynchronously. Wait for the asynchronous call to complete.
 e.wait();

}

//************************************ // Initialize the array from 0 to array_size - 1 //************************************ void InitializeArray(int *a, size_t size) {

 for (size_t i = 0; i < size; i++) a[i] = i;

}

//************************************ // Demonstrate vector add both in sequential on CPU and in parallel on device. //************************************ int main() {

 // Create device selector for the device of your interest.

1. if FPGA_EMULATOR

 // DPC++ extension: FPGA emulator selector on systems without FPGA card.
 INTEL::fpga_emulator_selector d_selector;

1. elif FPGA

 // DPC++ extension: FPGA selector on systems with FPGA card.
 INTEL::fpga_selector d_selector;

1. else

 // The default device selector will select the most performant device.
 default_selector d_selector;

1. endif

 try {
   queue q(d_selector, exception_handler);

   // Print out the device information used for the kernel code.
   std::cout << "Running on device: "
             << q.get_device().get_info<info::device::name>() << "\n";
   std::cout << "Vector size: " << array_size << "\n";

   // Create arrays with "array_size" to store input and output data. Allocate
   // unified shared memory so that both CPU and device can access them.
   int *a = malloc_shared<int>(array_size, q);
   int *b = malloc_shared<int>(array_size, q);
   int *sum_sequential = malloc_shared<int>(array_size, q);
   int *sum_parallel = malloc_shared<int>(array_size, q);

   if ((a == nullptr) || (b == nullptr) || (sum_sequential == nullptr) ||
       (sum_parallel == nullptr)) {
     if (a != nullptr) free(a, q);
     if (b != nullptr) free(b, q);
     if (sum_sequential != nullptr) free(sum_sequential, q);
     if (sum_parallel != nullptr) free(sum_parallel, q);

     std::cout << "Shared memory allocation failure.\n";
     return -1;
   }

   // Initialize input arrays with values from 0 to array_size - 1
   InitializeArray(a, array_size);
   InitializeArray(b, array_size);

   // Compute the sum of two arrays in sequential for validation.
   for (size_t i = 0; i < array_size; i++) sum_sequential[i] = a[i] + b[i];

   // Vector addition in DPC++.
   VectorAdd(q, a, b, sum_parallel, array_size);

   // Verify that the two arrays are equal.
   for (size_t i = 0; i < array_size; i++) {
     if (sum_parallel[i] != sum_sequential[i]) {
       std::cout << "Vector add failed on device.\n";
       return -1;
     }
   }

   int indices[]{0, 1, 2, (array_size - 1)};
   constexpr size_t indices_size = sizeof(indices) / sizeof(int);

   // Print out the result of vector add.
   for (int i = 0; i < indices_size; i++) {
     int j = indices[i];
     if (i == indices_size - 1) std::cout << "...\n";
     std::cout << "[" << j << "]: " << j << " + " << j << " = "
               << sum_sequential[j] << "\n";
   }

   free(a, q);
   free(b, q);
   free(sum_sequential, q);
   free(sum_parallel, q);
 } catch (exception const &e) {
   std::cout << "An exception is caught while adding two vectors.\n";
   std::terminate();
 }

 std::cout << "Vector add successfully completed on device.\n";
 return 0;

} ```

Math Mul

```//============================================================== // Copyright © 2020 Intel Corporation // // SPDX-License-Identifier: MIT // =============================================================

/**

* Matrix_mul multiplies two large matrices both the CPU and the offload device,
* then compares results. If the code executes on both CPU and the offload
* device, the name of the offload device and a success message are displayed.
*
* For comprehensive instructions regarding DPC++ Programming, go to
* https://software.intel.com/en-us/oneapi-programming-guide and search based on
* relevant terms noted in the comments.
*/

1. include <CL/sycl.hpp>
2. include <iostream>
3. include <limits>

// dpc_common.hpp can be found in the dev-utilities include folder. // e.g., $ONEAPI_ROOT/dev-utilities/<version>/include/dpc_common.hpp

1. include "dpc_common.hpp"

using namespace std; using namespace sycl;

/**

* Each element of the product matrix c[i][j] is computed from a unique row and
* column of the factor matrices, a[i][k] and b[k][j]
*/

// Matrix size constants. constexpr int m_size = 150 * 8; // Must be a multiple of 8. constexpr int M = m_size / 8; constexpr int N = m_size / 4; constexpr int P = m_size / 2;

/**

* Perform matrix multiplication on host to verify results from device.
*/

int VerifyResult(float (*c_back)[P]);

int main() {

 // Host memory buffer that device will write data back before destruction.
 float(*c_back)[P] = new float[M][P];

 // Intialize c_back
 for (int i = 0; i < M; i++)
   for (int j = 0; j < P; j++) c_back[i][j] = 0.0f;

 // Initialize the device queue with the default selector. The device queue is
 // used to enqueue kernels. It encapsulates all states needed for execution.
 try {
   queue q(default_selector{}, dpc_common::exception_handler);

   cout << "Device: " << q.get_device().get_info<info::device::name>() << "\n";

   // Create 2D buffers for matrices, buffer c is bound with host memory c_back

   buffer<float, 2> a_buf(range(M, N));
   buffer<float, 2> b_buf(range(N, P));
   buffer c_buf(reinterpret_cast<float *>(c_back), range(M, P));

   cout << "Problem size: c(" << M << "," << P << ") = a(" << M << "," << N
        << ") * b(" << N << "," << P << ")\n";

   // Using three command groups to illustrate execution order. The use of
   // first two command groups for initializing matrices is not the most
   // efficient way. It just demonstrates the implicit multiple command group
   // execution ordering.

   // Submit command group to queue to initialize matrix a
   q.submit([&](auto &h) {
     // Get write only access to the buffer on a device.
     accessor a(a_buf, h, write_only);

     // Execute kernel.
     h.parallel_for(range(M, N), [=](auto index) {
       // Each element of matrix a is 1.
       a[index] = 1.0f;
     });
   });

   // Submit command group to queue to initialize matrix b
   q.submit([&](auto &h) {
     // Get write only access to the buffer on a device
     accessor b(b_buf, h, write_only);

     // Execute kernel.
     h.parallel_for(range(N, P), [=](auto index) {
       // Each column of b is the sequence 1,2,...,N
       b[index] = index[0] + 1.0f;
     });
   });

   // Submit command group to queue to multiply matrices: c = a * b
   q.submit([&](auto &h) {
     // Read from a and b, write to c
     accessor a(a_buf, h, read_only);
     accessor b(b_buf, h, read_only);
     accessor c(c_buf, h, write_only);

     int width_a = a_buf.get_range()[1];

     // Execute kernel.
     h.parallel_for(range(M, P), [=](auto index) {
       // Get global position in Y direction.
       int row = index[0];
       // Get global position in X direction.
       int col = index[1];

       float sum = 0.0f;

       // Compute the result of one element of c
       for (int i = 0; i < width_a; i++) {
         sum += a[row][i] * b[i][col];
       }

       c[index] = sum;
     });
   });
 } catch (sycl::exception const &e) {
   cout << "An exception is caught while multiplying matrices.\n";
   terminate();
 }

 int result;
 cout << "Result of matrix multiplication using DPC++: ";
 result = VerifyResult(c_back);
 delete[] c_back;

 return result;

}

bool ValueSame(float a, float b) {

 return fabs(a - b) < numeric_limits<float>::epsilon();

}

int VerifyResult(float (*c_back)[P]) {

 // Check that the results are correct by comparing with host computing.
 int i, j, k;

 // 2D arrays on host side.
 float(*a_host)[N] = new float[M][N];
 float(*b_host)[P] = new float[N][P];
 float(*c_host)[P] = new float[M][P];

 // Each element of matrix a is 1.
 for (i = 0; i < M; i++)
   for (j = 0; j < N; j++) a_host[i][j] = 1.0f;

 // Each column of b_host is the sequence 1,2,...,N
 for (i = 0; i < N; i++)
   for (j = 0; j < P; j++) b_host[i][j] = i + 1.0f;

 // c_host is initialized to zero.
 for (i = 0; i < M; i++)
   for (j = 0; j < P; j++) c_host[i][j] = 0.0f;

 for (i = 0; i < M; i++) {
   for (k = 0; k < N; k++) {
     // Each element of the product is just the sum 1+2+...+n
     for (j = 0; j < P; j++) {
       c_host[i][j] += a_host[i][k] * b_host[k][j];
     }
   }
 }

 bool mismatch_found = false;

 // Compare host side results with the result buffer from device side: print
 // mismatched data 5 times only.
 int print_count = 0;

 for (i = 0; i < M; i++) {
   for (j = 0; j < P; j++) {
     if (!ValueSame(c_back[i][j], c_host[i][j])) {
       cout << "Fail - The result is incorrect for element: [" << i << ", "
            << j << "], expected: " << c_host[i][j]
            << ", but found: " << c_back[i][j] << "\n";
       mismatch_found = true;
       print_count++;
       if (print_count == 5) break;
     }
   }

   if (print_count == 5) break;
 }

 delete[] a_host;
 delete[] b_host;
 delete[] c_host;

 if (!mismatch_found) {
   cout << "Success - The results are correct!\n";
   return 0;
 } else {
   cout << "Fail - The results mismatch!\n";
   return -1;
 }

} ```