Computing 2D and 3D Trigonometric Transforms Using the 1D Transform in Intel® MKL

Submit New Article

January 18, 2011 10:30 AM PST


The current version of the Intel® Math Kernel Library (Intel® MKL) contains an example code for the use of the 1D Trigonometric Transformation (TT) function. This article shows how the use of that 1D TT function can be extended to perform 2D and 3D transforms.

Problem: Compute a 2D trigonometric transform using the 1D TT function.

To compute the 2D trigonometric transform of array f with N rows and N columns, we must compute the 1D trigonometric transform for each row and compute 1D trigonometric transform for each column.

The Intel MKL TT routines don’t support strided input data, so we must introduce data movement to our algorithm. The new algorithm:

Or, a parallel version of the algorithm for a 4-core processor could look like this:

Implementation note:

  • We recommend that you initialize, commit, and free Intel MKL TT functions outside of any parallel region.
  • You can use the same handle and dpar/spar for each thread in a parallel region.
  • You have to set all ipar_x/y[9] = "number of threads".
  • Also you can use the same ipar for each thread, but you will lose the warnings which are saved to ipar. 
d_init_trig_transform( &n_x, &tt_type_x, ipar_x, dpar_x, &ir );

ipar_x[9]= max_thread; // maximum number of threads used in calculations 
d_commit_trig_transform( f, &handle_x, ipar_x, dpar_x, &ir ); 
// Parallel for 
Loop over j {d_backward_trig_transform( f(:,j), &handle_x, ipar_x, dpar_x, &ir );} 
// end parallel 
free_trig_transform(&handle_x,ipar_x,&ir);

d_init_trig_transform( &n_y, &tt_type_y, ipar_y, dpar_y, &ir);
ipar_y[9]= max_thread; 
// ... 
d_commit_trig_transform( f, &handle_y, ipar_y, dpar_y, &ir );
// Parallel for 
Loop over i {d_backward_trig_transform( f(i,:), &handle_y, ipar_y, dpar_y, &ir );} 
// end parallel 
free_trig_transform( &handle_y, ipar_y, &ir);

Problem: Compute 3D trigonometric transform via 1D TT function.

To compute a 3D trigonometric transform of an array f with N*M*K elements we will need to compute the 1D trigonometric transform along each of the 3 dimensions.

Since the Intel MKL TT routines do not support strided input data we must use an algorithm with data movement:

 

A simple example is explained below.

Implementation note:

  • We recommend that you initialize, commit, and free Intel MKL TT functions outside of any parallel region.
  • You can use the same handle and dpar/spar for each thread in a parallel region.
  • You have to set all ipar_x/y/z[9] = "number of threads".
  • you can use the same ipar for each thread, but you will lose the warnings which are saved to ipar.


d_init_trig_transform( &n_x, &tt_type_x, ipar_x, dpar_x, &ir );

ipar_x[9]= max_thread; // maximum number of threads used in calculations 
d_commit_trig_transform( f, &handle_x, ipar_x, dpar_x, &ir ); 
// Parallel for 
Loop over j {d_backward_trig_transform( f(:,j,k), &handle_x, ipar_x, dpar_x,&ir);}
// end parallel 
free_trig_transform( &handle_x, ipar_x, &ir);
d_init_trig_transform(&n_y, &tt_type_y, ipar_y, dpar_y, &ir ); ipar_y[9]= max_thread; 
// ...

d_commit_trig_transform( f, &handle_y, ipar_y, dpar_y, &ir );

// Parallel for
Loop over i {d_backward_trig_transform( f(i,:, k), &handle_y, ipar_y, dpar_y,&ir);} 
// end parallel 
free_trig_transform( &handle_y, ipar_y, &ir); 

d_init_trig_transform( &n_z, &tt_type_z, ipar_z, dpar_z,&ir ); 
ipar_z[9]= max_thread; 
// ... 
d_commit_trig_transform( f, &handle_z, ipar_z, dpar_z, &ir); 
// Parallel for 
Loop over i { d_backward_trig_transform( f(i,j,:), &handle_z, ipar_z, dpar_z, &ir);} 
// end parallel 
free_trig_transform( &handle_z, ipar_z, &ir);


Do you need more help?