DAPPLE Example: dirichlet

Source code:

// @TITLE "dirichlet.cc -- parallel Dirichlet program" // // dirichlet -- a parallel Dirichlet relaxation program // written using the Data-Parallel Programming Library for Education (DAPPLE) // // The Dirichlet function is really just a 4-nearest-neighbors // averaging stencil, applied repeatedly over many iterations // to a matrix whose boundary conditions are fixed; eventually, // the values converge to a steady state. In this example we // pick some arbitrary boundary values so that two opposite corners // are high, the other two corners are low, and the other boundary // values are linearly interpolated between the corners. // // Note that we precompute a boolean vector that remembers which // points are internal, to simplify the loop. Also, note the use of // a MAX reduction to determine convergence. // // David Kotz 1994 // $Id: dirichlet.cc,v 1.5 94/09/29 12:18:02 dfk Exp Locker: dfk $ #include #include #include #include "dapple.h" // we'll work with NxN matrices const int N = 5; const float epsilon = 0.05; // which row or column am I? useful for VP-specific math const intMatrix VProw(N, N, IdentityR); const intMatrix VPcol(N, N, IdentityC); int main(int argc, char **argv) { floatMatrix M(N,N); // current values of matrix floatMatrix Mprev(N,N); // previous values of matrix booleanMatrix isInternal(N,N); // false only on the edges of the matrix int iteration = 0; // initialize the matrix // internal points should be zero, but boundary points could be anything // I chose something simple. Note the corners are initialized with the // top and bottom rows, although the function I've chosen for the boundary // values is consistent at the corners. ifp (VProw == 0) { // top boundary isInternal = false; M = VPcol + 1; // 1..N } else ifp (VProw == N-1) { // bottom boundary isInternal = false; M = N - VPcol; // N..1 } else ifp (VPcol == 0) { // left boundary isInternal = false; M = VProw + 1; // 1..N } else ifp (VPcol == N-1) { // right boundary isInternal = false; M = N - VProw; // N..1 } else { // internal points isInternal = true; M = 0; } cout << "initial state:" << endl; cout << M << endl; // iterate until the maximum change is epsilon or less do { Mprev = M; // remember previous value ifp (isInternal) { M = ( shift(M, 0, -1) // right neighbor + shift(M, 0, 1) // left neighbor + shift(M, -1, 0) // top neighbor + shift(M, 1, 0)); // bottom neighbor M = M / 4; } cout << "after iteration " << ++iteration << ":" << endl; cout << M << endl; } while (max_value(M - Mprev) > epsilon); return(0); }

Demonstration:

dirichlet initial state: 1 2 3 4 5 2 0 0 0 4 3 0 0 0 3 4 0 0 0 2 5 4 3 2 1 after iteration 1: 1 2 3 4 5 2 1 0.75 2 4 3 0.75 0 0.75 3 4 2 0.75 1 2 5 4 3 2 1 after iteration 2: 1 2 3 4 5 2 1.375 1.5 2.375 4 3 1.5 0.75 1.5 3 4 2.375 1.5 1.375 2 5 4 3 2 1 after iteration 3: 1 2 3 4 5 2 1.75 1.875 2.75 4 3 1.875 1.5 1.875 3 4 2.75 1.875 1.75 2 5 4 3 2 1 after iteration 4: 1 2 3 4 5 2 1.9375 2.25 2.9375 4 3 2.25 1.875 2.25 3 4 2.9375 2.25 1.9375 2 5 4 3 2 1 after iteration 5: 1 2 3 4 5 2 2.125 2.4375 3.125 4 3 2.4375 2.25 2.4375 3 4 3.125 2.4375 2.125 2 5 4 3 2 1 after iteration 6: 1 2 3 4 5 2 2.21875 2.625 3.21875 4 3 2.625 2.4375 2.625 3 4 3.21875 2.625 2.21875 2 5 4 3 2 1 after iteration 7: 1 2 3 4 5 2 2.3125 2.71875 3.3125 4 3 2.71875 2.625 2.71875 3 4 3.3125 2.71875 2.3125 2 5 4 3 2 1 after iteration 8: 1 2 3 4 5 2 2.35938 2.8125 3.35938 4 3 2.8125 2.71875 2.8125 3 4 3.35938 2.8125 2.35938 2 5 4 3 2 1 after iteration 9: 1 2 3 4 5 2 2.40625 2.85938 3.40625 4 3 2.85938 2.8125 2.85938 3 4 3.40625 2.85938 2.40625 2 5 4 3 2 1 after iteration 10: 1 2 3 4 5 2 2.42969 2.90625 3.42969 4 3 2.90625 2.85938 2.90625 3 4 3.42969 2.90625 2.42969 2 5 4 3 2 1