Now, take a look at three different integration methods: naïve summation, two-point trapezoid, and a three-point Simpson's rule. The following is the expression for the example:
1 #include <blitz/vector.h>
2 usingnamespace blitz;
3
4 int main()
5 {
6 // Integrate the expression
7
8 // to estimate the error function erf(x) on the interval [0,1].
9 //
10 // Three methods are compared:
11 // 1. Naive summation
12 // 2. Two-point trapezoid
13 // 3. Three-point Simpson's
14
15 constint numSamples = 1024; // Number of samples
16 double dt = 1.0 / numSamples; // Distance between samples
17 Range R(0, numSamples - 1); // Index set for sample// vectors
18 Range displayRange(0, 900, 100); // Indices to output to// screen
19
20 // For comparison purposes, use the built-in erf(x) function.
21 Vector<double> z = erf(R * dt);
22 cout << " Built-in: " << z(displayRange+1) << endl;
23
24
25 // Naive summation
26 Vector<double> z1 = accumulate(M_2_SQRTPI *
exp(-sqr(R * dt)) * dt);
27 cout << " Naive: " << z1(displayRange+1) << endl;
28
29 // Calculate the rms error
30 double error1 = sqrt(mean(sqr(z1 - z)));
31 cout << "RMS Error: " << error1 << endl;
32
33
34 // For the following rules, it's easier to work with a
35 // sampled integrand.
36 Vector<double> a = M_2_SQRTPI * exp(-sqr(R * dt));
37
38
39 // Two-point trapezoid
40 Range J(1, numSamples-1);
41 Vector<double> z2 = accumulate(dt / 2.0 * (a(J) + a(J-1)));
42
43 cout << " Trapezoid: " << z2(displayRange) << endl;
44 cout << "RMS Error: " << sqrt(mean(sqr(z2-z(J)))) << endl;
45
46
47 // Three-point Simpson's (parabolic)
48 Range I(1, numSamples-2);
49 Vector<double> z3 = accumulate(dt / 6.0 * (a(I-1) + 4 *
a(I) + a(I+1)));
50
51 cout << " Parabolic: " << z3(displayRange) << endl;
52 cout << " RMS Error: " << sqrt(mean(sqr(z3 - z(I)))) << endl;
53
54
55 return 0;
56 }
57
Lines 26, 41, and 49 use the accumulate functor. A functor describes an action that should be taken on each member of a list. The action usually will vary according to the list member's type.
Parallel Computation
With multi-core CPUs on the horizon, everyone is thinking about how to exploit the extra cores by adding computation threads. Although Blitz++ can be used for parallel computing, it was not designed primarily for that. For this reason, you may want to investigate some other available libraries such as POOMA before choosing to implement a parallel code using Blitz++. Using threads on Blitz++ may require a library recompile, depending on platform.
In threadsafe mode, Blitz++ array reference counts are safeguarded by a mutex. By default, pthread mutexes are used. Blitz++ does not do locking for every array element access because it would result in an unacceptable performance hit. Responsibility for appropriate synchronization is squarely on your shoulders.
Tough Numeric Problems? No Problem
You've just barely scratched the surface of what Blitz++ can do for you. Some of the topics this article couldn't cover include sophisticated array expressions, stencils, slicing, and indirection. Give it a try because Blitz++ is ready, willing, and able to handle your toughest numeric problems.
The Numerical Recipes series is a classic among numerical algorithm books. Thoroughly self-contained, it proceeds from mathematical and theoretical considerations to actual, practical computer routines. In addition to C++, Numerical Recipes books are also available for C, Fortran 77, Fortran 90, and Pascal.
About the Author
Victor Volkman has been writing for C/C++ Users Journal and other programming journals since the late 1980s. He is a graduate of Michigan Tech and a faculty advisor board member for Washtenaw Community College CIS department. Volkman is the editor of numerous books, including C/C++ Treasure Chest and is the owner of Loving Healing Press. He can help you in your quest for open source tools and libraries, just drop an e-mail to sysop@HAL9K.com.
About the Author
Victor Volkman has been writing for C/C++ Users Journal and other programming journals since the late 1980s. He is a graduate of Michigan Tech and a faculty advisor board member for Washtenaw Community College CIS department. Volkman is the editor of numerous books, including C/C++ Treasure Chest and is the owner of Loving Healing Press. He can help you in your quest for open source tools and libraries, just drop an e-mail to sysop@HAL9K.com.
Add www.codeguru.com to your favorites Add www.codeguru.com to your browser search box IE 7 | Firefox 2.0 | Firefox 1.5.xReceive news via our XML/RSS feed
IE 7
Firefox 2.0
