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