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