Using Powerful C++ Extension to Solve Scientific Problems Share your comment!

The Cilk Plus extensions include array notation, which goes beyond the C++ standard, to help you easily write vectorized code. Among its many uses, the powerful mechanism can be a very valuable tool for scientific programmers and engineers. Jeff Cogswell explains how. 

In my previous blog posting, I introduced the concept of the array notation that’s unique to the Intel compiler. As with the three Cilk Plus keywords, array notations are additions to the compiler that go beyond the C++ standard. They make use of vectorization and multi-core parallelism to provide for highly parallel and optimized code.

After spending a lot of time with scientific programmers and engineers, I’ve realized that many of them could really benefit from using the array notation in their scientific modeling, but actually weren’t aware of it. So, in this blog and the next, I’m going to tackle a few scientific applications.

Linear Algebra

At the heart of physics is linear algebra. Many scientific problems can be modeled through the concepts of vectors and matrices. Elementary physics courses start out with teaching how vectors can work together.

In the most basic form, a vector is simply a set of numbers in a specific order. In physics, a vector can represent a force, including the direction and how strong it is (called the magnitude). When you push on something, you provide force that will cause the object to move in a certain direction. The force also has a strength to it. You can show the force graphically with an arrow pointing in a certain direction, with the arrow’s length representing how strong the force is. If you start the vector at the origin of a three-dimensional space, the end of the vector will have a three-point coordinate. Those three numbers can be used to fully represent the vector. Thus, a force might have a vector with these numbers:

{ 1.0, 2.0, 3.0 }

When forces operate in conjunction, you can add their components to get the resulting force. So if you have another force that looks like this:

{ 7.0, 9.0, 11.0 }

And these two forces are together pushing an object, the final resulting force is simply the addition of the components:

{ 1.0 + 7.0, 2.0 + 9.0, 3.0 + 11.0 }

Using the array notation I discussed in the previous blog, we can easily add these up and get a new vector in a single operation. Here’s some code that does exactly that:

float A[3] = { 1.0, 2.0, 3.0 };
float B[3] = { 7.0, 9.0, 11.0 };
float C[3];
C[:] = A[:] + B[:];

printf("{ %.2f, %.2f, %.2f }", C[0], C[1], C[2]);

The fourth line is key. This line adds all the elements of the arrays, using vectorization when possible. (We also have to make sure our data is aligned, which I’ll cover in the next blog. Note that I’m also using the handy printf statement so I can easily format my output, instead of stringing together a bunch of cout insertions.)

Dot Products

Another common problem in engineering is the concept of a dot product. With this we can lead to more advanced matrix operations, which are also on the plate for next time. But meanwhile, a dot product is simply an operation where you take the corresponding components of two vectors and multiply them together to get a new vector. Then take that vector and sum the components. For example, if our two vectors are as above, we first multiply the corresponding terms like so:

{ 1.0 x 7.0, 2.0 x 9.0, 3.0 x 11.0 }

to get a vector like this:

{ 7.0, 18.0, 33.0 }

And then we add them up:

7.0 + 18.0 + 33.0 = 58.0

Doing this with Cilk Plus is almost deceptively easy. The Cilk Plus runtime library includes a function called __sec_reduce_add, which uses reducers to add up numbers in parallel. We can pass into it the result of an array operation, like so:

float x = __sec_reduce_add(A[:] * B[:]);

where A and B are as defined in the code up above. Notice I’m doing it all in a single step. The operation

A[:] * B[:]

creates a vector whose components are the products of the corresponding components in A and B. Then I pass this vector to __sec_reduc_add, which adds up the members of the vector. And just like that, we have our total. Here’s some code that does it:

float A[3] = { 1.0, 2.0, 3.0 };
float B[3] = { 7.0, 9.0, 11.0 };
float x = __sec_reduce_add(A[:] * B[:]);
printf("%.2f", x);


The Array Notation is very powerful and certainly useful in scientific and engineering applications. Next time, we’ll look at how to align the data for maximum performance, as well as more advanced matrix operations. Meanwhile, share your thoughts in the comments below regarding your own scientific applications and how you’ve found Cilk Plus to help.

Posted on July 15, 2013 by Jeff Cogswell, Slashdot Media Contributing Editor