DAPPLE Example: bisection
Source code:
// @TITLE "bisection.cc -- parallel bisection"
//
// bisection -- a parallel bisection program
// written using the Data-Parallel Programming Library for Education (DAPPLE)
//
// This program finds roots of an equation using a parallel bisection method.
// In the sequential version of this method, we are given two points
// and the assumption that a root does exist between
// those two points. We narrow the range on each iteration by testing the
// midpoint and comparing the sign of the equation at the endpoints and
// at the midpoint, much like a binary search.
// In this parallel version of the algorithm, the range is broken into
// N-1 intervals (for N processors) in each iteration; the processors
// then compare their value with the one in the processor to their right,
// to see which pair of proessors bracket the root. They redefine the
// range for the next iteration.
// Note that the function f(x) is defined as a vector function so that
// many values can be computed (for different x values) in parallel,
// with Y = f(X). Otherwise, with a plain f(x) function, it could have
// been written Y = apply(f, X).
//
// David Kotz 1994
// $Id: bisection.cc,v 2.13 94/09/29 12:17:55 dfk Exp Locker: dfk $
#include
#include
#include
#include "dapple.h"
// #define DEBUG // turns on some debugging output
const int N = 11;
const float a0 = -5;
const float b0 = 50;
const float epsilon_x = 0.0001; // stop when (b-a) < epsilon_x
// which VP am I? useful for VP-specific math
const intVector VP(N, Identity);
inline floatVector f(const floatVector& x) { return x*x - 10; }
int
main(int argc, char **argv)
{
float a = a0, b = b0;
float interval_size;
floatVector X(N), Y(N);
while ((b-a) >= epsilon_x) {
// divide it into N-1 intervals
interval_size = (b-a) / (N-1);
cout << "a=" << a << "; b=" << b <<
"; interval=" << interval_size << endl;
X = a + interval_size * floatVector(VP);
Y = f(X);
#ifdef DEBUG
cout << X << endl;
cout << Y << endl;
cout << endl;
#endif
if (any(Y == 0)) { // did we find the actual root?
ifp (Y == 0) // found a root, where is it?
// there should only be one VP active here
a = b = first(X);
} else { // no root found; bracket the root
ifp (VP > 0 && Y * shift(Y, 1) < 0) {
// there should only be one VP active here
a = first(shift(X, 1));
b = first(X);
}
}
}
cout << "The root is between " << a << " and " << b << endl;
return(0);
}
Demonstration:
bisection
a=-5; b=50; interval=5.5
a=-5; b=0.5; interval=0.55
a=-3.35; b=-2.8; interval=0.055
a=-3.185; b=-3.13; interval=0.00550001
a=-3.163; b=-3.1575; interval=0.000550008
a=-3.16245; b=-3.1619; interval=5.50032e-05
The root is between -3.16228 and -3.16223